comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Various improvements of the helium code (mostly by lwalewsk)

Browse source 
repaired and adjusted to comply with current cp2k by (fuhl,cschran)


svn-origin-rev: 16409
This commit is contained in:
Harald Forbert 2015-12-07 09:46:51 +00:00
parent 887e240406
commit ded76cc2c1
21 changed files with 7877 additions and 4435 deletions
Download patch file Download diff file Expand all files Collapse all files

29
src/common/bibliography.F
View file

@ -74,7 +74,9 @@ MODULE bibliography
Khaliullin2007, Khaliullin2008, Merlot2014, Lin2009, Lin2013, &
QUIP_ref, DelBen2015, Souza2002, Umari2002, Stengel2009, &
Luber2014, Berghold2011, DelBen2015b, Campana2009, &
Schiffmann2015, Bruck2014, Rappe1992, Ceriotti2012, Ceriotti2010
Schiffmann2015, Bruck2014, Rappe1992, Ceriotti2012, &
Ceriotti2010, Walewski2014
CONTAINS
@ -2378,6 +2380,31 @@ CONTAINS
"ER"),&
DOI="10.1103/PhysRevB.87.041108")
CALL add_reference(key=Walewski2014, ISI_record=s2a(&
"PT J",&
"AU Walewski, L",&
" Forbert, H",&
" Marx, D",&
"AF Walewski, Lukasz",&
" Forbert, Harald",&
" Marx, Dominik",&
"TI Reactive path integral quantum simulations of molecules solvated in",&
" superfluid helium",&
"SO COMPUTER PHYSICS COMMUNICATIONS",&
"DE Cryochemistry; Superfluidity; Helium; Solvation; Path integrals",&
"NR 153",&
"SN 0010-4655",&
"PD MAR",&
"PY 2014",&
"VL 185",&
"IS 3",&
"BP 884",&
"EP 899",&
"DI 10.1016/j.cpc.2013.12.011",&
"UT WOS:000331919100022",&
"ER"),&
DOI="10.1016/j.cpc.2013.12.011")
CALL add_reference(key=Delben2013, ISI_record=s2a(&
"PT J",&
"AU Del Ben, M",&

1581
src/helium_common.F
View file

File diff suppressed because it is too large Load diff

1200
src/helium_io.F
View file

File diff suppressed because it is too large Load diff

187
src/helium_types.F
View file

@ -1,187 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2015 CP2K developers group !
!-----------------------------------------------------------------------------!
! *****************************************************************************
!> \brief Data types representing superfluid helium
!> \author hforbert
!> \date 2009-01-01
!> \par History
!> extracted helium_solvent_type from pint_types.F [lwalewski]
! *****************************************************************************
MODULE helium_types
USE cell_types, ONLY: cell_type
USE cp_log_handling, ONLY: cp_logger_type
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: dp,&
int_8
USE parallel_rng_types, ONLY: rng_stream_type
USE splines_types, ONLY: spline_data_p_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
REAL(kind=dp), PARAMETER, PUBLIC :: he_mass = 4.00263037059764_dp ! 4He mass in [u]
! Energy contributions - symbolic names for indexing energy arrays
INTEGER, PARAMETER, PUBLIC :: e_id_total = 1, &
e_id_potential = 2, &
e_id_kinetic = 3, &
e_id_interact = 4, &
e_id_thermo = 5, &
e_id_virial = 6
! Number of energy contributions for static array allocation
INTEGER, PARAMETER, PUBLIC :: e_num_ids = 10
PUBLIC :: helium_solvent_type
! *****************************************************************************
!> \brief data of solvent helium
!> \note
!> ref_count - reference count of this data structure
!> id_nr - identification number of this data structure
!> input - input data structure (the whole tree)
!> num_env - number of independent He environments
!> num_env_restart - number of He environments present in the restart file
!> periodic - true if bulk liquid helium in periodic box
!> cell_shape - unit cell shape for PBC calculations
!> cell_size - size of the periodic box for helium only
!> cell_size_inv - 1 / cell_size (inverse)
!> cell_m_inv - invrse of the unit cell vectors' matrix
!> tau - 1/(k_B T p) with T being the He temperature, p - number of beads
!> density - helium density for free bulk in box
!> e_corr - potential correction energy due to finite box
!> hb2m - hbar squared over m for helium
!> atoms - number of helium atoms
!> beads - number of helium beads per atom (needs to be an integer
!> multiple of the solvated molecule bead number)
!> pdx - pair density expansion max exponent
!> maxcycle - maximum cyclic permutation change to attempt
!> m_value - cycle length sampled with different probability than other lengths
!> m_ratio - probability ratio betw m_value and other possible values of m
!> bisection - power of 2 number for bisection algorithm
!> bisctlog2 - log2(bisection)
!> relrot - relative rotation in imaginary time with
!> respect to the normal system/starting configuration
!> iter_norot - number of iterations to try for a given
!> imaginary time slice rotation
!> iter_rot - number of rotations to try
!> total number of iterations is iter_norot*iter_rot
!> iter - number of iterations tried so far
!> first_step - first step, restarted from MOTION%PINT%ITERATION (def val =0)
!> current_step - first_step + number of steps performed so far
!> accepts - number of accepted new configurations
!> num_accepted - average number of accepted permutations of a given length
!> on a given Levy level, plus one additional level which
!> counts # of trials, REAL(BISCTLOG2+2, MAX_PERM_CYCLE)
!> num_accepted(1,l) - # of trials for perm length l
!> num_accepted(2,l) - # of selected perms of length l
!> num_accepted(3,l) - # of perms of length l accepted at level 1
!> average over He environments/processors
!> wnumber_inst - winding number (instantaneous) REAL(3)
!> wnumber_avrg - winding number (averaged) REAL(3)
!> sdensity_inst - superfluid density (instantaneous) REAL
!> sdensity_avrg - superfluid density (averaged) REAL
!> pos - position of the helium atoms REAL(3,ATOMS,BEADS)
!> work - same dimensions as pos
!> tmatrix - ? permutation probability related
!> pmatrix - ? permutation probability related
!> [use might change/new ones added/etc]
!> pweight - ? permutation probability related
!> ptable - proposed cyclic permutation
!> (dimension max_cycle)
!> permutation - current permutation state INTEGER(ATOMS)
!> iperm - inverse of the current permutation state INTEGER(ATOMS)
!> plength_avrg - permutation length probability distribution REAL(ATOMS)
!> plength_inst - instantaneous permutation length probability REAL(ATOMS)
!> atom_plength - length of the permutation cycle the atom belongs to INTEGER(ATOMS)
!> uij - pair density matrix coefficients (action)
!> eij - pair density matrix coefficients (energy)
!> bead_ratio - ratio of helium beads to system beads
!> rng_stream_uniform - random numbers from uniform distribution
!> rng_stream_gaussian - random numbers from gaussian distribution
!> rdf_delr - delta r for RDF
!> rdf_maxr - maximum r for RDF
!> solute_present- switch the interactions with the solute on or off
!> solute_atoms - number of solute atoms (=pint_env%ndim/3)
!> solute_beads - number of solute beads (=pint_env%p)
!> solute_element- element names of solute atoms (NDIM/3)
!> solute_number - number of solute atoms of different atomic kinds
!> solute_index - indices of solute atoms sorted by atomic kinds
!> solute_cell - dimensions of the solvated system cell (a,b,c)
!> (not needed now and should be removed at some point)
!> force_avrg - averaged forces exerted by He solvent on the solute (P,NDIM)
!> force_inst - instantaneous forces exerted by He on the solute (P, NDIM)
!> General purpose temporary arrays
!> (performance measure: no allocation in frequently called subroutines)
!> rtmp_3_np_1d - real, 1D temp array (3 * NUM_ENV)
!> rtmp_p_ndim_1d - real, 1D temp array (P * NDIM)
!> rtmp_p_ndim_np_1d - real, 1D temp array (P * NDIM * NUM_ENV)
!> rtmp_p_ndim_2d - real, 2D temp array (P, NDIM)
!> rtmp_3_atoms_beads_1d - real, 1D temp array (3 * ATOMS * BEADS)
!> ltmp_3_atoms_beads_3d - logical, 3D temp array (3,ATOMS,BEADS) - same as pos
!> itmp_atoms_1d - integer, 1D temp array (ATOMS) - same as permutation
!> itmp_atoms_np_1d - integer, 1D temp array (ATOMS*NUM_ENV)
!> ltmp_atoms_1d - logical, 1D temp array (ATOMS) - for unpacking permutation
!> \author hforbert
! *****************************************************************************
TYPE helium_solvent_type
INTEGER :: id_nr,ref_count
TYPE(section_vals_type), POINTER :: input
TYPE(cp_logger_type), POINTER :: logger
LOGICAL :: periodic
INTEGER :: cell_shape, num_env, num_env_restart, m_value
REAL(kind=dp) :: cell_size, cell_size_inv, m_ratio
REAL(kind=dp), DIMENSION(3,3) :: cell_m_inv
REAL(kind=dp) :: tau, density, e_corr, hb2m, pweight
INTEGER :: atoms, beads, pdx, maxcycle, relrot, iter_norot, iter_rot
INTEGER :: bisection, bisctlog2
INTEGER :: bead_ratio, num_steps, first_step, last_step, current_step
INTEGER(kind=int_8) :: accepts
INTEGER, DIMENSION(:), POINTER :: ptable
INTEGER, DIMENSION(:), POINTER :: permutation
INTEGER, DIMENSION(:), POINTER :: iperm
REAL(kind=dp), DIMENSION(:,:,:), POINTER :: pos
REAL(kind=dp), DIMENSION(:,:,:), POINTER :: work
REAL(kind=dp), DIMENSION(:,:), POINTER :: tmatrix
REAL(kind=dp), DIMENSION(:,:), POINTER :: pmatrix
REAL(kind=dp), DIMENSION(:,:), POINTER :: ipmatrix
INTEGER, DIMENSION(:,:), POINTER :: nmatrix
TYPE (spline_data_p_type), DIMENSION(:,:), POINTER :: uij
TYPE (spline_data_p_type), DIMENSION(:,:), POINTER :: eij
REAL(kind=dp), DIMENSION(e_num_ids) :: energy_inst, energy_avrg
REAL(kind=dp), DIMENSION(3) :: wnumber_inst, wnumber_avrg, origin
REAL(kind=dp) :: sdensity_inst, sdensity_avrg
REAL(kind=dp), DIMENSION(:,:), POINTER :: num_accepted
REAL(kind=dp), DIMENSION(:), POINTER :: plength_avrg, plength_inst
INTEGER, DIMENSION(:), POINTER :: atom_plength
INTEGER :: rdf_nbin
REAL(kind=dp) :: rdf_delr, rdf_maxr
REAL(kind=dp), DIMENSION(:), POINTER :: rdf_inst, rdf_avrg
INTEGER :: rho_nbin, rho_minb, rho_iweight, rho_num
REAL(kind=dp) :: rho_delr, rho_maxr
REAL(kind=dp), DIMENSION(:,:,:,:), POINTER :: rho_inst, rho_avrg, rho_rstr
LOGICAL :: rho_restart, rho_present
TYPE(rng_stream_type), POINTER :: rng_stream_gaussian, rng_stream_uniform
! variables that describe the solvated molecular system
LOGICAL :: solute_present
INTEGER :: solute_atoms, solute_beads
CHARACTER(LEN=2), DIMENSION(:), POINTER :: solute_element
INTEGER, DIMENSION(:), POINTER :: solute_number
INTEGER, DIMENSION(:,:), POINTER :: solute_index
TYPE(cell_type), POINTER :: solute_cell
REAL(KIND=dp), DIMENSION(:,:), POINTER :: force_avrg, force_inst
INTEGER, DIMENSION(:), POINTER :: itmp_atoms_1d, itmp_atoms_np_1d
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_3_np_1d, rtmp_p_ndim_1d, rtmp_p_ndim_np_1d
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_3_atoms_beads_1d, rtmp_3_atoms_beads_np_1d
REAL(KIND=dp), DIMENSION(:,:), POINTER :: rtmp_p_ndim_2d
LOGICAL, DIMENSION(:,:,:), POINTER :: ltmp_3_atoms_beads_3d
LOGICAL, DIMENSION(:), POINTER :: ltmp_atoms_1d
END TYPE helium_solvent_type
END MODULE helium_types

35
src/input_constants.F
View file

@ -288,6 +288,41 @@ MODULE input_constants
helium_cell_shape_cube=1,&
helium_cell_shape_octahedron=2
! helium-solute interaction potentials
INTEGER, PARAMETER, PUBLIC :: helium_solute_intpot_none=0,&
helium_solute_intpot_mwater=1
! helium force selection
INTEGER, PARAMETER, PUBLIC :: helium_forces_average = 0, &
helium_forces_last = 1
! superfluid density, estimators
INTEGER, PARAMETER, PUBLIC :: estimator_none=1,&
estimator_parea=2,&
estimator_wnumber=3, &
estimator_pcycle=4, &
estimator_weighted=5, &
estimator_normalization=6
! superfluid density, denominator selector
INTEGER, PARAMETER, PUBLIC :: denominator_unity=1,&
denominator_natoms=2,&
denominator_inertia=3,&
denominator_rperp2=4
! helium-sampling algorithms
INTEGER, PARAMETER, PUBLIC :: helium_sampling_ceperley=0
! distribution types for sampling path lengths
INTEGER, PARAMETER, PUBLIC :: helium_mdist_singlev=1,&
helium_mdist_uniform=2,&
helium_mdist_linear=3,&
helium_mdist_quadratic=4,&
helium_mdist_exponential=5,&
helium_mdist_gaussian=6
INTEGER, PARAMETER, PUBLIC :: perm_plain=1,&
perm_cycle=2
! Free Energy methods
INTEGER, PARAMETER, PUBLIC :: do_fe_meta=0,&
do_fe_ui=1,&

2620
src/motion/helium_common.F Normal file
View file

File diff suppressed because it is too large Load diff

721
src/motion/helium_interactions.F
View file

@ -12,13 +12,12 @@ MODULE helium_interactions
USE helium_common, ONLY: helium_eval_expansion,&
helium_pbc
USE helium_types, ONLY: e_id_interact,&
e_id_kinetic,&
e_id_potential,&
e_id_thermo,&
e_id_total,&
e_id_virial,&
helium_solvent_type
USE helium_types, ONLY: &
e_id_interact, e_id_kinetic, e_id_potential, e_id_thermo, e_id_total, &
e_id_virial, helium_solvent_type, hid_chlorine, hid_hydrogen, &
hid_oxygen, int_arr_ptr
USE input_constants, ONLY: helium_solute_intpot_mwater,&
helium_solute_intpot_none
USE kinds, ONLY: dp
USE physcon, ONLY: angstrom,&
kelvin
@ -35,7 +34,8 @@ MODULE helium_interactions
PUBLIC :: helium_calc_energy
PUBLIC :: helium_solute_e_f
PUBLIC :: helium_bead_solute_e
PUBLIC :: helium_bead_solute_e_f
PUBLIC :: helium_intpot_scan
CONTAINS
@ -161,21 +161,18 @@ vkin = 0.0_dp
RETURN
END SUBROUTINE helium_calc_energy
! ***************************************************************************
!> \brief Calculate helium-solute interaction energy and forces contribution
!> from one helium bead of one helium particle and the corresponding
!> solute time slice.
!> \param pint_env ...
!> \brief Calculate general helium-solute interaction energy (and forces)
!> between one helium bead and the corresponding solute time slice.
!> \param pint_env path integral environment
!> \param helium ...
!> \param helium_part_index ...
!> \param helium_slice_index ...
!> \param helium_r_opt ...
!> \param energy ...
!> \param force ...
!> \date 2009-12-11
!> \param helium_part_index helium particle index
!> \param helium_slice_index helium time slice index
!> \param helium_r_opt explicit helium bead coordinates (optional)
!> \param energy calculated energy
!> \param force calculated force (if requested)
!> \author Lukasz Walewski
!> \note This is a stub, Lennard-Jones potential with no warranty. It should
!> be replaced with realistic potential for real-world simulations.
! *****************************************************************************
SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
helium_slice_index, helium_r_opt, energy, force)
@ -188,24 +185,14 @@ vkin = 0.0_dp
INTENT(IN), OPTIONAL :: helium_r_opt
REAL(KIND=dp), INTENT(OUT) :: energy
REAL(KIND=dp), DIMENSION(:, :), &
INTENT(OUT) :: force
INTENT(OUT), OPTIONAL, POINTER :: force
CHARACTER(len=*), PARAMETER :: routineN = 'helium_bead_solute_e_f', &
routineP = moduleN//':'//routineN
INTEGER :: hbeads, hi, i, ig, &
num_chlorine, num_hydrogen, &
num_oxygen, qi
REAL(KIND=dp) :: d, d2, dd, ep, eps, s1, s2, &
sig
REAL(KIND=dp), DIMENSION(3) :: dr, helium_r, solute_r
CPASSERT(ASSOCIATED(helium))
CPASSERT(ASSOCIATED(pint_env))
num_chlorine = helium%solute_number(1)
num_oxygen = helium%solute_number(2)
num_hydrogen = helium%solute_number(3)
INTEGER :: hbeads, hi, qi
REAL(KIND=dp), DIMENSION(3) :: helium_r
REAL(KIND=dp), DIMENSION(:), POINTER :: my_force
hbeads = helium%beads
! helium bead index that is invariant wrt the rotations
@ -220,78 +207,47 @@ vkin = 0.0_dp
helium_r(:) = helium%pos(:,helium_part_index,helium_slice_index)
END IF
energy = 0.0_dp
force(:,:) = 0.0_dp
SELECT CASE (helium%solute_interaction)
sig = 2.69_dp ! 1.4 Angstrom
eps = 60.61e-6_dp ! 19 K
s1 = 0.0_dp
DO i = 1, num_hydrogen
ig = helium%solute_index(3,i)-1 ! global hydrogen index (3 == H)
solute_r(1) = pint_env%x(qi,3*ig+1)
solute_r(2) = pint_env%x(qi,3*ig+2)
solute_r(3) = pint_env%x(qi,3*ig+3)
dr(:) = solute_r(:) - helium_r(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
s2 = 24.0_dp*eps*dd*(2.0_dp*dd-1.0_dp)/d2
force(qi,3*ig+1) = force(qi,3*ig+1) + s2*dr(1)
force(qi,3*ig+2) = force(qi,3*ig+2) + s2*dr(2)
force(qi,3*ig+3) = force(qi,3*ig+3) + s2*dr(3)
END DO ! i = 1, num_hydrogen
energy = energy + s1
sig = 5.01_dp ! 2.6 Angstrom
eps = 79.17e-6_dp ! 25 K
s1 = 0.0_dp
DO i = 1, num_chlorine
ig = helium%solute_index(1,i)-1 ! global chlorine index (1 == Cl)
solute_r(1) = pint_env%x(qi,3*ig+1)
solute_r(2) = pint_env%x(qi,3*ig+2)
solute_r(3) = pint_env%x(qi,3*ig+3)
dr(:) = solute_r(:) - helium_r(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
s2 = 24.0_dp*eps*dd*(2.0_dp*dd-1.0_dp)/d2
force(qi,3*ig+1) = force(qi,3*ig+1) + s2*dr(1)
force(qi,3*ig+2) = force(qi,3*ig+2) + s2*dr(2)
force(qi,3*ig+3) = force(qi,3*ig+3) + s2*dr(3)
END DO ! i = 1, num_chlorine
energy = energy + s1
CASE (helium_solute_intpot_mwater)
IF (PRESENT(force)) THEN
force(:,:) = 0.0_dp
my_force => force(qi,:)
CALL helium_intpot_model_water( &
pint_env%x(qi,:), &
helium%solute_i, &
helium, &
helium_r,&
energy, &
my_force &
)
ELSE
CALL helium_intpot_model_water( &
pint_env%x(qi,:), &
helium%solute_i, &
helium, &
helium_r,&
energy &
)
END IF
sig = 5.01_dp ! 2.6 Angstrom
eps = 104.5e-6_dp ! 33 K
s1 = 0.0_dp
DO i = 1, num_oxygen
ig = helium%solute_index(2,i)-1 ! global oxygen index (2 == O)
solute_r(1) = pint_env%x(qi,3*ig+1)
solute_r(2) = pint_env%x(qi,3*ig+2)
solute_r(3) = pint_env%x(qi,3*ig+3)
dr(:) = solute_r(:) - helium_r(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
s2 = 24.0_dp*eps*dd*(2.0_dp*dd-1.0_dp)/d2
force(qi,3*ig+1) = force(qi,3*ig+1) + s2*dr(1)
force(qi,3*ig+2) = force(qi,3*ig+2) + s2*dr(2)
force(qi,3*ig+3) = force(qi,3*ig+3) + s2*dr(3)
END DO ! i = 1, num_chlorine
energy = energy + s1
CASE (helium_solute_intpot_none)
energy = 0.0_dp
IF (PRESENT(force)) THEN
force(:,:) = 0.0_dp
END IF
CASE DEFAULT
END SELECT
RETURN
END SUBROUTINE helium_bead_solute_e_f
! ***************************************************************************
!> \brief Calculate total helium-solute interaction energy and forces.
!> \param pint_env - path integral environment
@ -308,7 +264,7 @@ vkin = 0.0_dp
CHARACTER(len=*), PARAMETER :: routineN = 'helium_solute_e_f', &
routineP = moduleN//':'//routineN
INTEGER :: i, ibead, ipart, j
INTEGER :: ia, ib, jb, jc
REAL(KIND=dp) :: my_energy
REAL(KIND=dp), DIMENSION(:, :), POINTER :: force
@ -320,14 +276,14 @@ vkin = 0.0_dp
! calculate the total interaction energy and gradients between the
! solute and the helium, sum over all beads of all He particles
DO ipart = 1, helium%atoms
DO ibead = 1, helium%beads
CALL helium_bead_solute_e_f(pint_env, helium, ipart, &
ibead, energy=my_energy, force=helium%rtmp_p_ndim_2d)
energy = energy + my_energy
DO i = 1, pint_env%p
DO j = 1, pint_env%ndim
force(i,j) = force(i,j) + helium%rtmp_p_ndim_2d(i,j)
DO ia = 1, helium%atoms
DO ib = 1, helium%beads
CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, &
energy=my_energy, force=helium%rtmp_p_ndim_2d)
energy = energy + my_energy
DO jb = 1, pint_env%p
DO jc = 1, pint_env%ndim
force(jb,jc) = force(jb,jc) + helium%rtmp_p_ndim_2d(jb,jc)
END DO
END DO
END DO
@ -336,113 +292,6 @@ vkin = 0.0_dp
RETURN
END SUBROUTINE helium_solute_e_f
! ***************************************************************************
!> \brief Calculate helium-solute interaction energy contribution from one
!> helium bead of one helium particle and the corresponding solute
!> time slice.
!> \param pint_env ...
!> \param helium ...
!> \param helium_part_index ...
!> \param helium_slice_index ...
!> \param helium_r_opt ...
!> \param energy ...
!> \date 2009-12-11
!> \author Lukasz Walewski
!> \note This is a stub, Lennard-Jones potential with no warranty. It should
!> be replaced with realistic potential for real-world simulations.
! *****************************************************************************
SUBROUTINE helium_bead_solute_e(pint_env, helium, helium_part_index, &
helium_slice_index, helium_r_opt, energy)
TYPE(pint_env_type), POINTER :: pint_env
TYPE(helium_solvent_type), POINTER :: helium
INTEGER, INTENT(IN) :: helium_part_index, &
helium_slice_index
REAL(KIND=dp), DIMENSION(3), &
INTENT(IN), OPTIONAL :: helium_r_opt
REAL(KIND=dp), INTENT(OUT) :: energy
INTEGER :: hbeads, hi, i, ig, &
num_chlorine, num_hydrogen, &
num_oxygen, qi
REAL(KIND=dp) :: d, d2, dd, ep, eps, s1, sig
REAL(KIND=dp), DIMENSION(3) :: dr, helium_r, solute_r
num_chlorine = helium%solute_number(1)
num_oxygen = helium%solute_number(2)
num_hydrogen = helium%solute_number(3)
hbeads = helium%beads
! helium bead index that is invariant wrt the rotations
hi = MOD(helium_slice_index-1+hbeads+helium%relrot,hbeads) + 1
! solute bead index that belongs to hi helium index
qi = ((hi-1)*pint_env%p)/hbeads+1
! coordinates of the helium bead
IF (PRESENT(helium_r_opt)) THEN
helium_r(:) = helium_r_opt(:)
ELSE
helium_r(:) = helium%pos(:,helium_part_index,helium_slice_index)
END IF
energy = 0.0_dp
sig = 2.69_dp ! 1.4 Angstrom
eps = 60.61e-6_dp ! 19 K
s1 = 0.0_dp
DO i = 1, num_hydrogen
ig = helium%solute_index(3,i)-1 ! global hydrogen index (3 == H)
solute_r(1) = pint_env%x(qi,3*ig+1)
solute_r(2) = pint_env%x(qi,3*ig+2)
solute_r(3) = pint_env%x(qi,3*ig+3)
dr(:) = solute_r(:) - helium_r(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1)+dr(2)*dr(2)+dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
END DO ! i = 1, num_hydrogen
energy = energy + s1
sig = 5.01_dp ! 2.6 Angstrom
eps = 79.17e-6_dp ! 25 K
s1 = 0.0_dp
DO i = 1, num_chlorine
ig = helium%solute_index(1,i)-1 ! global chlorine index (1 == Cl)
solute_r(1) = pint_env%x(qi,3*ig+1)
solute_r(2) = pint_env%x(qi,3*ig+2)
solute_r(3) = pint_env%x(qi,3*ig+3)
dr(:) = solute_r(:) - helium_r(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1)+dr(2)*dr(2)+dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
END DO ! i = 1, num_chlorine
energy = energy + s1
sig = 5.01_dp ! 2.6 Angstrom
eps = 104.5e-6_dp ! 33 K
s1 = 0.0_dp
DO i = 1, num_oxygen
ig = helium%solute_index(2,i)-1 ! global oxygen index (2 == O)
solute_r(1) = pint_env%x(qi,3*ig+1)
solute_r(2) = pint_env%x(qi,3*ig+2)
solute_r(3) = pint_env%x(qi,3*ig+3)
dr(:) = solute_r(:) - helium_r(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1)+dr(2)*dr(2)+dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
END DO ! i = 1, num_oxygen
energy = energy + s1
RETURN
END SUBROUTINE helium_bead_solute_e
! ***************************************************************************
!> \brief Calculate total helium-solute interaction energy.
@ -457,15 +306,14 @@ vkin = 0.0_dp
TYPE(helium_solvent_type), POINTER :: helium
REAL(KIND=dp), INTENT(OUT) :: energy
INTEGER :: ibead, ipart
INTEGER :: ia, ib
REAL(KIND=dp) :: my_energy
energy = 0.0_dp
DO ipart = 1, helium%atoms
DO ibead = 1, helium%beads
CALL helium_bead_solute_e(pint_env, helium, ipart, &
ibead, energy=my_energy)
DO ia = 1, helium%atoms
DO ib = 1, helium%beads
CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, energy=my_energy)
energy = energy + my_energy
END DO
END DO
@ -473,6 +321,214 @@ vkin = 0.0_dp
RETURN
END SUBROUTINE helium_solute_e
! ***************************************************************************
!> \brief Calculate l-th Legendre polynomial P_{l}(x) at point x
!> \param x ...
!> \param n ...
!> \retval Pl ...
! *****************************************************************************
FUNCTION Pl(x,n)
REAL(KIND=dp), INTENT(IN) :: x
INTEGER, INTENT(IN) :: n
REAL(KIND=dp) :: Pl
INTEGER :: k
REAL(KIND=dp) :: pln(0:n)
pln(0) = 1.0_dp
pln(1) = x
IF (n <= 1) THEN
Pl = pln(n)
ELSE
DO k=1,n-1
pln(k+1) = ((2.0*k+1.0)*x*pln(k) - REAL(k,dp)*pln(k-1))/(REAL(k+1,dp))
END DO
Pl = pln(n)
END IF
RETURN
END FUNCTION Pl
! ***************************************************************************
!> \brief Scan the helium-solute interaction energy within the periodic cell
!> \param pint_env ...
!> \param helium ...
!> \date 2014-01-22
!> \author Lukasz Walewski
! *****************************************************************************
SUBROUTINE helium_intpot_scan(pint_env, helium)
TYPE(pint_env_type), POINTER :: pint_env
TYPE(helium_solvent_type), POINTER :: helium
CHARACTER(len=*), PARAMETER :: routineN = 'helium_intpot_scan', &
routineP = moduleN//':'//routineN
INTEGER :: handle, ic, ix, iy, iz, nbin
LOGICAL :: wrapped
REAL(KIND=dp) :: delr, my_en, ox, oy, oz
REAL(kind=dp), DIMENSION(3) :: pbc1, pbc2, pos
CALL timeset(routineN,handle)
nbin = helium%rho_nbin
delr = helium%rho_delr
ox = helium%center(1) - helium%rho_maxr / 2.0_dp
oy = helium%center(2) - helium%rho_maxr / 2.0_dp
oz = helium%center(3) - helium%rho_maxr / 2.0_dp
DO ix = 1, nbin
DO iy = 1, nbin
DO iz = 1, nbin
! put the probe in the center of the current voxel
pos(:) = (/ox+(ix-0.5_dp)*delr, oy+(iy-0.5_dp)*delr, oz+(iz-0.5_dp)*delr/)
! calc interaction energy for the current probe position
helium%pos(:,1,1) = pos(:)
CALL helium_bead_solute_e_f(pint_env,helium,1,1,energy=my_en)
! check if the probe fits within the unit cell
pbc1(:) = pos(:) - helium%center
pbc2(:) = pbc1(:)
CALL helium_pbc( helium, pbc2 )
wrapped = .FALSE.
DO ic = 1, 3
IF (ABS(pbc1(ic)-pbc2(ic)) .GT. 10.0_dp * EPSILON(0.0_dp)) THEN
wrapped = .TRUE.
END IF
END DO
! set the interaction energy value
IF (wrapped) THEN
helium%rho_inst(1,ix,iy,iz) = 0.0_dp
ELSE
helium%rho_inst(1,ix,iy,iz) = my_en
END IF
END DO
END DO
END DO
CALL timestop(handle)
RETURN
END SUBROUTINE helium_intpot_scan
! ***************************************************************************
!> \brief Calculate model helium-solute interaction energy and forces
!> between one helium bead and the corresponding solute time
!> slice asuming water solute.
!> \param solute_x solute positions ARR(3*NATOMS)
!> \param solute_i solute indices, mapping from instances of a given element
!> to global atom indices
!> \param helium only needed for helium_pbc call at the moment
!> \param helium_x helium bead position ARR(3)
!> \param energy calculated interaction energy
!> \param force ...
!> \author Felix Uhl
! *****************************************************************************
SUBROUTINE helium_intpot_model_water(solute_x, solute_i, helium, helium_x, energy, force)
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: solute_x
TYPE(int_arr_ptr), DIMENSION(:), &
INTENT(IN) :: solute_i
TYPE(helium_solvent_type), POINTER :: helium
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: helium_x
REAL(KIND=dp), INTENT(OUT) :: energy
REAL(KIND=dp), DIMENSION(:), &
INTENT(OUT), OPTIONAL, POINTER :: force
CHARACTER(LEN=*), PARAMETER :: routineN = 'helium_intpot_model_water', &
routineP = moduleN//':'//routineN
INTEGER :: i, ig, num_chlorine, &
num_hydrogen, num_oxygen
REAL(KIND=dp) :: d, d2, dd, ep, eps, s1, s2, &
sig
REAL(KIND=dp), DIMENSION(3) :: dr, solute_r
IF (ASSOCIATED(solute_i(hid_chlorine)%iap)) THEN
num_chlorine = SIZE(solute_i(hid_chlorine)%iap)
ELSE
num_chlorine = 0
END IF
IF (ASSOCIATED(solute_i(hid_oxygen)%iap)) THEN
num_oxygen = SIZE(solute_i(hid_oxygen)%iap)
ELSE
num_oxygen = 0
END IF
IF (ASSOCIATED(solute_i(hid_hydrogen)%iap)) THEN
num_hydrogen = SIZE(solute_i(hid_hydrogen)%iap)
ELSE
num_hydrogen = 0
END IF
energy = 0.0_dp
IF (PRESENT(force)) THEN
force(:) = 0.0_dp
END IF
sig = 2.69_dp ! 1.4 Angstrom
eps = 60.61e-6_dp ! 19 K
s1 = 0.0_dp
DO i = 1, num_hydrogen
ig = solute_i(hid_hydrogen)%iap(i)-1 ! global hydrogen index (3 == H)
solute_r(1) = solute_x(3*ig+1)
solute_r(2) = solute_x(3*ig+2)
solute_r(3) = solute_x(3*ig+3)
dr(:) = solute_r(:) - helium_x(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
s2 = 24.0_dp*eps*dd*(2.0_dp*dd-1.0_dp)/d2
IF (PRESENT(force)) THEN
force(3*ig+1) = force(3*ig+1) + s2*dr(1)
force(3*ig+2) = force(3*ig+2) + s2*dr(2)
force(3*ig+3) = force(3*ig+3) + s2*dr(3)
END IF
END DO ! i = 1, num_hydrogen
energy = energy + s1
sig = 5.01_dp ! 2.6 Angstrom
eps = 104.5e-6_dp ! 33 K
s1 = 0.0_dp
DO i = 1, num_oxygen
ig = solute_i(hid_oxygen)%iap(i)-1 ! global oxygen index (2 == O)
solute_r(1) = solute_x(3*ig+1)
solute_r(2) = solute_x(3*ig+2)
solute_r(3) = solute_x(3*ig+3)
dr(:) = solute_r(:) - helium_x(:)
CALL helium_pbc( helium, dr )
d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
d = SQRT(d2)
dd = (sig/d)**6
ep = 4.0_dp*eps*dd*(dd-1.0_dp)
s1 = s1 + ep
s2 = 24.0_dp*eps*dd*(2.0_dp*dd-1.0_dp)/d2
IF (PRESENT(force)) THEN
force(3*ig+1) = force(3*ig+1) + s2*dr(1)
force(3*ig+2) = force(3*ig+2) + s2*dr(2)
force(3*ig+3) = force(3*ig+3) + s2*dr(3)
END IF
END DO ! i = 1, num_chlorine
energy = energy + s1
RETURN
END SUBROUTINE helium_intpot_model_water
! ***************************************************************************
!> \brief Helium-helium pair interaction potential.
!> \param r ...
@ -498,6 +554,11 @@ vkin = 0.0_dp
RETURN
END FUNCTION helium_vij
#if 0
! this block is currently turned off
! ***************************************************************************
!> \brief Helium-helium pair interaction potential's derivative.
!> \param r ...
@ -530,4 +591,224 @@ vkin = 0.0_dp
RETURN
END FUNCTION helium_d_vij
! *****************************************************************************
!> \brief ...
!> \param helium ...
!> \param n ...
!> \param i ...
!> \retval res ...
! *****************************************************************************
FUNCTION helium_atom_action(helium,n,i) RESULT(res)
TYPE(helium_solvent_type), POINTER :: helium
INTEGER, INTENT(IN) :: n, i
REAL(KIND=dp) :: res
INTEGER :: c, j
REAL(KIND=dp) :: r(3), rp(3), s, t
s = 0.0_dp
t = 0.0_dp
IF (n < helium%beads) THEN
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,i,n+1)
END DO
CALL helium_pbc( helium, r )
t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
DO j = 1, i-1
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,j,n)
rp(c) = helium%pos(c,i,n+1) - helium%pos(c,j,n+1)
END DO
CALL helium_pbc( helium, r )
CALL helium_pbc( helium, rp )
s = s + helium_eval_expansion(helium,r,rp,helium%uij,1)
END DO
DO j = i+1, helium%atoms
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,j,n)
rp(c) = helium%pos(c,i,n+1) - helium%pos(c,j,n+1)
END DO
CALL helium_pbc( helium, r )
CALL helium_pbc( helium, rp )
s = s + helium_eval_expansion(helium,r,rp,helium%uij,1)
END DO
ELSE
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,helium%permutation(i),1)
END DO
CALL helium_pbc( helium, r )
t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
DO j = 1, i-1
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,j,n)
rp(c) = helium%pos(c,helium%permutation(i),1) - helium%pos(c,helium%permutation(j),1)
END DO
CALL helium_pbc( helium, r )
CALL helium_pbc( helium, rp )
s = s + helium_eval_expansion(helium,r,rp,helium%uij,1)
END DO
DO j = i+1, helium%atoms
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,j,n)
rp(c) = helium%pos(c,helium%permutation(i),1) - helium%pos(c,helium%permutation(j),1)
END DO
CALL helium_pbc( helium, r )
CALL helium_pbc( helium, rp )
s = s + helium_eval_expansion(helium,r,rp,helium%uij,1)
END DO
END IF
t = t / (2.0_dp*helium%tau*helium%hb2m)
s = s * 0.5_dp
res = s+t
RETURN
END FUNCTION helium_atom_action
! *****************************************************************************
!> \brief ...
!> \param helium ...
!> \param n ...
!> \retval res ...
! *****************************************************************************
FUNCTION helium_link_action(helium,n) RESULT(res)
TYPE(helium_solvent_type), POINTER :: helium
INTEGER, INTENT(IN) :: n
REAL(KIND=dp) :: res
INTEGER :: c, i, j
REAL(KIND=dp) :: r(3), rp(3), s, t
s = 0.0_dp
t = 0.0_dp
IF (n < helium%beads) THEN
DO i = 1, helium%atoms
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,i,n+1)
END DO
CALL helium_pbc( helium, r )
t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
DO j = 1, i-1
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,j,n)
rp(c) = helium%pos(c,i,n+1) - helium%pos(c,j,n+1)
END DO
CALL helium_pbc( helium, r )
CALL helium_pbc( helium, rp )
s = s + helium_eval_expansion(helium,r,rp,helium%uij,1)
END DO
END DO
ELSE
DO i = 1, helium%atoms
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,helium%permutation(i),1)
END DO
CALL helium_pbc( helium, r )
t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
DO j = 1, i-1
DO c = 1, 3
r(c) = helium%pos(c,i,n) - helium%pos(c,j,n)
rp(c) = helium%pos(c,helium%permutation(i),1) - helium%pos(c,helium%permutation(j),1)
END DO
CALL helium_pbc( helium, r )
CALL helium_pbc( helium, rp )
s = s + helium_eval_expansion(helium,r,rp,helium%uij,1)
END DO
END DO
END IF
t = t / (2.0_dp*helium%tau*helium%hb2m)
res = s+t
RETURN
END FUNCTION helium_link_action
! *****************************************************************************
!> \brief ...
!> \param helium ...
!> \retval res ...
! *****************************************************************************
FUNCTION helium_total_action(helium) RESULT(res)
TYPE(helium_solvent_type), POINTER :: helium
REAL(KIND=dp) :: res
INTEGER :: i
REAL(KIND=dp) :: s
s = 0.0_dp
DO i = 1, helium%beads
s = s + helium_link_action(helium,i)
END DO
res = s
RETURN
END FUNCTION helium_total_action
! *****************************************************************************
!> \brief ...
!> \param helium ...
!> \param part ...
!> \param ref_bead ...
!> \param delta_bead ...
!> \param d ...
! *****************************************************************************
SUBROUTINE helium_delta_pos(helium,part,ref_bead,delta_bead,d)
TYPE(helium_solvent_type), POINTER :: helium
INTEGER, INTENT(IN) :: part, ref_bead, delta_bead
REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: d
INTEGER :: b, bead, db, nbead, np, p
REAL(KIND=dp), DIMENSION(3) :: r
b = helium%beads
d(:) = 0.0_dp
IF (delta_bead > 0) THEN
bead = ref_bead
p = part
db = delta_bead
DO
IF (db < 1) EXIT
nbead = bead + 1
np = p
IF (nbead > b) THEN
nbead = nbead - b
np = helium%permutation(np)
END IF
r(:) = helium%pos(:,p,bead) - helium%pos(:,np,nbead)
CALL helium_pbc( helium, r )
d(:) = d(:) + r(:)
bead = nbead
p = np
db = db-1
END DO
ELSEIF ( delta_bead < 0) THEN
bead = ref_bead
p = part
db = delta_bead
DO
IF (db >= 0) EXIT
nbead = bead - 1
np = p
IF (nbead < 1) THEN
nbead = nbead + b
np = helium%iperm(np)
END IF
r(:) = helium%pos(:,p,bead) - helium%pos(:,np,nbead)
CALL helium_pbc( helium, r )
d(:) = d(:) + r(:)
bead = nbead
p = np
db = db + 1
END DO
END IF
RETURN
END SUBROUTINE helium_delta_pos
#endif
END MODULE helium_interactions

1675
src/motion/helium_io.F Normal file
View file

File diff suppressed because it is too large Load diff

1862
src/motion/helium_methods.F
View file

File diff suppressed because it is too large Load diff

1194
src/motion/helium_sampling.F
View file

File diff suppressed because it is too large Load diff

398
src/motion/helium_types.F Normal file
View file

@ -0,0 +1,398 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2015 CP2K developers group !
!-----------------------------------------------------------------------------!
! *****************************************************************************
!> \brief Data types representing superfluid helium
!> \author hforbert
!> \date 2009-01-01
!> \par History
!> extracted helium_solvent_type from pint_types.F [lwalewski]
! *****************************************************************************
MODULE helium_types
USE cell_types, ONLY: cell_type
USE cp_log_handling, ONLY: cp_logger_type
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: default_string_length,&
dp,&
int_8
USE parallel_rng_types, ONLY: rng_stream_type
USE splines_types, ONLY: spline_data_p_type
#include "../base/base_uses.f90"
IMPLICIT NONE
PRIVATE
LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_types'
REAL(KIND=dp), PARAMETER, PUBLIC :: he_mass = 4.00263037059764_dp !< 4He mass in [u]
!> Energy contributions - symbolic names for indexing energy arrays
INTEGER, PARAMETER, PUBLIC :: &
e_id_total = 1, &
e_id_potential = 2, &
e_id_kinetic = 3, &
e_id_interact = 4, &
e_id_thermo = 5, &
e_id_virial = 6
!> Number of energy contributions for static array allocation
INTEGER, PARAMETER, PUBLIC :: e_num_ids = 10
!> number of density function identifiers
INTEGER, PARAMETER, PUBLIC :: rho_num = 5
!> density function identifier names
INTEGER, PARAMETER, PUBLIC :: &
rho_atom_number = 1, &
rho_projected_area = 2,&
rho_winding_number = 3,&
rho_winding_cycle = 4,&
rho_moment_of_inertia = 5
!> number of atom type identifiers
INTEGER, PARAMETER, PUBLIC :: hid_num = 4
!> atom type identifier names (for helium-solute interactions)
INTEGER, PARAMETER, PUBLIC :: &
hid_chlorine = 1, &
hid_oxygen = 2,&
hid_hydrogen = 3, &
hid_carbon = 4
!> derived data types
PUBLIC :: helium_solvent_type
PUBLIC :: int_arr_ptr
PUBLIC :: real_arr_ptr
!> functions
PUBLIC :: helium_destroy_int_arr_ptr
PUBLIC :: helium_destroy_real_arr_ptr
! ***************************************************************************
!> \brief Vector type useful for averaging
!> \author Lukasz Walewski
!> \date 2014-09-09
! ***************************************************************************
TYPE helium_vector_type
!> instantaneous value
REAL(KIND=dp), DIMENSION(3) :: inst
!> accumulated value
REAL(KIND=dp), DIMENSION(3) :: accu
!> running average
REAL(KIND=dp), DIMENSION(3) :: ravr
!> restarted value
REAL(KIND=dp), DIMENSION(3) :: rstr
END TYPE helium_vector_type
! ***************************************************************************
!> \brief data structure for solvent helium
!> \author hforbert
! ***************************************************************************
TYPE helium_solvent_type
INTEGER :: id_nr !< identification number of this data structure
INTEGER :: ref_count !< reference count of this data structure
TYPE(section_vals_type), POINTER :: input !< input data structure (the whole tree)
TYPE(cp_logger_type), POINTER :: logger
INTEGER :: num_env !< number of He environments in runtime
INTEGER :: num_env_restart !< number of He environments in restart file
INTEGER :: atoms !< number of atoms
INTEGER :: beads !< number of beads per atom (needs to be an integer multiple of the solute's number of beads)
INTEGER :: bead_ratio !< ratio of helium beads to system beads
REAL(KIND=dp) :: density !< helium density for free bulk in box
! some useful constants
!
REAL(KIND=dp) :: hb2m !< hbar squared over m for 4He in CP2k units
REAL(KIND=dp) :: tau !< 1/(k_B T p) with T - He temperature, p - number of beads
REAL(KIND=dp) :: wpref !< prefactor for calculating superfluid fraction from <(M*W)^2>
REAL(KIND=dp) :: apref !< prefactor for calculating superfluid fraction from <A^2/I_c>
! PBC related
!
LOGICAL :: periodic !< true if bulk liquid helium in periodic box
INTEGER :: cell_shape !< unit cell shape for PBC calculations
REAL(KIND=dp) :: cell_size !< size of the periodic box (helium only)
REAL(KIND=dp) :: cell_size_inv !< 1/cell_size (inverse)
REAL(KIND=dp), DIMENSION(3,3) :: cell_m !< the unit cell vectors' matrix
REAL(KIND=dp), DIMENSION(3,3) :: cell_m_inv !< invrse of the unit cell vectors' matrix
REAL(KIND=dp), DIMENSION(3) :: origin !< origin of the cell (first voxel position)
REAL(KIND=dp) :: droplet_radius !< radius of the droplet
REAL(KIND=dp), DIMENSION(3) :: center !< COM of solute (if present) or center of periodic cell (if periodic) or COM of helium
! sampling parameters
INTEGER :: sampling_method
INTEGER :: iter_norot !< number of iterations to try for a given imaginary time slice rotation (num inner MC loop iters)
INTEGER :: iter_rot !< number of rotations to try (total number of iterations is iter_norot*iter_rot) (num outer MC loop iters)
!
INTEGER :: maxcycle !< maximum cyclic permutation change to attempt
INTEGER :: m_dist_type !< distribution from which the cycle length m is sampled
INTEGER :: m_value !< cycle length sampled with different probability than other lengths
REAL(KIND=dp) :: m_ratio !< probability ratio betw m_value and other possible values of m
!
INTEGER :: relrot !< relative rotation in imaginary time wrt normal system/starting configuration
INTEGER :: bisection !< power of 2 number for bisection algorithm
INTEGER :: bisctlog2 !< log2(bisection)
REAL(KIND=dp) :: e_corr !< potential correction energy due to finite box
INTEGER :: pdx !< pair density expansion max exponent
! MC step counters
!
INTEGER :: num_steps !< number of iterations in the current run
INTEGER :: first_step !< first step, restarted from MOTION%PINT%ITERATION (default value: 0)
INTEGER :: last_step
INTEGER :: current_step !< first_step + number of steps performed so far
! helium variables
!
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: pos !< position of the helium atoms DIM(3,atoms,beads)
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: work !< same dimensions as pos
!
INTEGER, DIMENSION(:), POINTER :: permutation !< current permutation state DIM(atoms)
INTEGER, DIMENSION(:), POINTER :: iperm !< inverse of the current permutation state DIM(atoms)
INTEGER, DIMENSION(:), POINTER :: ptable !< proposed cyclic permutation, DIM(max_cycle)
INTEGER(KIND=int_8) :: accepts !< number of accepted new configurations
!
REAL(KIND=dp), DIMENSION(:,:), POINTER :: tmatrix !< ? permutation probability related
REAL(KIND=dp), DIMENSION(:,:), POINTER :: pmatrix !< ? permutation probability related [use might change/new ones added/etc]
REAL(KIND=dp) :: pweight !< ? permutation probability related
REAL(KIND=dp), DIMENSION(:,:), POINTER :: ipmatrix
INTEGER, DIMENSION(:,:), POINTER :: nmatrix
TYPE (spline_data_p_type), DIMENSION(:,:), POINTER :: uij !< pair density matrix coefficients (action)
TYPE (spline_data_p_type), DIMENSION(:,:), POINTER :: eij !< pair density matrix coefficients (energy)
! calculated properties
!
REAL(KIND=dp), DIMENSION(e_num_ids) :: energy_inst !< energy contributions (instantaneous)
REAL(KIND=dp), DIMENSION(e_num_ids) :: energy_avrg !< energy contributions (averaged)
TYPE(helium_vector_type) :: wnumber !< winding number
TYPE(helium_vector_type) :: wnmber2 !< winding number squared
TYPE(helium_vector_type) :: proarea !< projected area
TYPE(helium_vector_type) :: prarea2 !< projected area squared
TYPE(helium_vector_type) :: mominer !< moment of inertia
INTEGER :: averages_iweight !< weight for restarted averages
LOGICAL :: averages_restarted !< flag indicating whether the averages have been restarted
!
INTEGER :: rdf_nbin !< number of bins for RDF
INTEGER :: rdf_iweight !< weight for restarted RDF
INTEGER :: rho_iweight !< weight for restarted RHO
INTEGER :: rdf_num !< number of RDF estimators
INTEGER :: rdf_num_ctr !< number of centers for RDF calc
REAL(KIND=dp) :: rdf_delr !< delta r for RDF
REAL(KIND=dp) :: rdf_maxr !< maximum r for RDF
REAL(KIND=dp), DIMENSION(:), POINTER :: rdf_centers !< positions of RDF centers
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: rdf_inst !< RDF (instantaneous/tmp array)
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: rdf_rstr !< RDF (restarted)
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: rdf_accu !< RDF (accumulated for one run)
LOGICAL :: rdf_present
!
INTEGER :: rho_nbin
INTEGER :: rho_num_act !< actual number of density estimators
INTEGER :: rho_num_min_len_wdg !< number of optional estimators based on winding cycles
INTEGER :: rho_num_min_len_non!< number of optional estimators based on non-winding cycles
INTEGER :: rho_num_min_len_all!< number of optional estimators based on all cycles
INTEGER, DIMENSION(:), POINTER :: rho_min_len_wdg_vals !< minimum lengths of winding cycles
INTEGER, DIMENSION(:), POINTER :: rho_min_len_non_vals !< minimum lengths of non-winding cycles
INTEGER, DIMENSION(:), POINTER :: rho_min_len_all_vals !< minimum lengths of all cycles
REAL(KIND=dp) :: rho_delr, rho_maxr
REAL(KIND=dp), DIMENSION(:,:,:,:), POINTER :: rho_inst
REAL(KIND=dp), DIMENSION(:,:,:,:), POINTER :: rho_rstr
REAL(KIND=dp), DIMENSION(:,:,:,:), POINTER :: rho_accu
LOGICAL :: rho_present
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: rho_incr !< increment for density bining
TYPE(density_properties_type), DIMENSION(:), POINTER :: rho_property
REAL(KIND=dp), DIMENSION(:,:), POINTER :: num_accepted !< average number of accepted permutations of a given length
!! on a given Levy level, plus one additional level which
!! counts # of trials, REAL(BISCTLOG2+2, MAX_PERM_CYCLE)
!! num_accepted(1,l) - # of trials for perm length l
!! num_accepted(2,l) - # of selected perms of length l
!! num_accepted(3,l) - # of perms of length l accepted at level 1
!! average over He environments/processors
REAL(KIND=dp), DIMENSION(:), POINTER :: plength_avrg !< permutation length probability distribution DIM(atoms)
REAL(KIND=dp), DIMENSION(:), POINTER :: plength_inst !< instantaneous permutation length probability DIM(atoms)
INTEGER, DIMENSION(:), POINTER :: atom_plength !< length of the permutation cycle the atom belongs to DIM(atoms)
TYPE(rng_stream_type), POINTER :: rng_stream_uniform !< random number stream with uniform distribution
TYPE(rng_stream_type), POINTER :: rng_stream_gaussian !< random number stream with gaussian distribution
! variables related to solvated molecular system
!
LOGICAL :: solute_present !< switch the interactions with the solute on or off
INTEGER :: solute_atoms !< number of solute atoms (pint_env%ndim/3)
INTEGER :: solute_beads !< number of solute beads (pint_env%p)
INTEGER :: get_helium_forces !< parameter to determine whether the average or last MC force should be taken to MD
CHARACTER(LEN=2), DIMENSION(:), POINTER :: solute_element !< element names of solute atoms (pint_env%ndim/3)
TYPE(int_arr_ptr), DIMENSION(:), POINTER :: solute_i !< indices of solute atoms sorted by atomic kinds
TYPE(cell_type), POINTER :: solute_cell !< dimensions of the solvated system cell (a,b,c) (not needed now and should be removed at some point)
REAL(KIND=dp), DIMENSION(:,:), POINTER :: force_avrg !< averaged forces exerted by He solvent on the solute DIM(p,ndim)
REAL(KIND=dp), DIMENSION(:,:), POINTER :: force_inst !< instantaneous forces exerted by He on the solute (p,ndim)
CHARACTER(LEN=2), DIMENSION(:), POINTER :: ename
INTEGER, DIMENSION(:), POINTER :: eid
INTEGER :: enum
INTEGER :: solute_interaction
LOGICAL :: interaction_pot_scan !< wheather to perform solute-helium interaction scan
! temporary arrays for optimization
!
INTEGER, DIMENSION(:), POINTER :: itmp_atoms_1d !< DIM(atoms) - same as permutation
INTEGER, DIMENSION(:), POINTER :: itmp_atoms_np_1d !< DIM(atoms*num_env)
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_3_np_1d !< DIM(3*num_env)
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_p_ndim_1d !< DIM(p*ndim)
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_p_ndim_np_1d !< DIM(p*ndim*num_env)
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_3_atoms_beads_1d !< DIM(3*atoms*beads)
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_3_atoms_beads_np_1d
REAL(KIND=dp), DIMENSION(:,:), POINTER :: rtmp_p_ndim_2d !< DIM(p,ndim)
LOGICAL, DIMENSION(:,:,:), POINTER :: ltmp_3_atoms_beads_3d !< DIM(3,atoms,beads) - same as pos
LOGICAL, DIMENSION(:), POINTER :: ltmp_atoms_1d !< DIM(atoms) - for unpacking the permutation
END TYPE helium_solvent_type
! ***************************************************************************
!> \brief Container type for properties of a helium density function
!> \author Lukasz Walewski
!> \date 2014-09-09
! ***************************************************************************
TYPE density_properties_type
!> name of this density function
CHARACTER(len=default_string_length) :: name
!> flag indicating whether this function should be calculated
LOGICAL :: is_calculated
!> number of components that this function is composed of
INTEGER :: num_components
!> suffixes for the filenames storing components of this function
CHARACTER(len=default_string_length), DIMENSION(:), POINTER :: filename_suffix
!> component names
CHARACTER(len=default_string_length), DIMENSION(:), POINTER :: component_name
!> indices locating the components of this function in the global density arrays
INTEGER, DIMENSION(:), POINTER :: component_index
END TYPE density_properties_type
! ***************************************************************************
!> \brief A pointer to an integer array, data type to be used in arrays of
!> pointers.
!> \author Lukasz Walewski
!> \date 2013-12-11
! ***************************************************************************
TYPE int_arr_ptr
INTEGER, DIMENSION(:), POINTER :: iap
END TYPE int_arr_ptr
! ***************************************************************************
!> \brief A pointer to a real array, data type to be used in arrays of
!> pointers.
!> \author Lukasz Walewski
!> \date 2013-12-11
! ***************************************************************************
TYPE real_arr_ptr
REAL(KIND=dp), DIMENSION(:), POINTER :: rap
END TYPE real_arr_ptr
CONTAINS
! ***************************************************************************
!> \brief Deallocate all arrays pointed to by the pointers stored in the
!> integer pointer array
!> \param int_arr_p ...
!> \date 2013-12-12
!> \author Lukasz Walewski
! *****************************************************************************
SUBROUTINE helium_destroy_int_arr_ptr(int_arr_p)
TYPE(int_arr_ptr), DIMENSION(:), POINTER :: int_arr_p
CHARACTER(len=*), PARAMETER :: routineN = 'helium_destroy_int_arr_ptr', &
routineP = moduleN//':'//routineN
INTEGER :: ip
! deallocate memory used by each component of the pointer array
DO ip = 1, SIZE(int_arr_p)
IF (ASSOCIATED(int_arr_p(ip)%iap)) THEN
DEALLOCATE(int_arr_p(ip)%iap)
END IF
END DO
! deallocate the memory used for pointer array
IF (ASSOCIATED(int_arr_p)) THEN
DEALLOCATE(int_arr_p)
END IF
RETURN
END SUBROUTINE helium_destroy_int_arr_ptr
! ***************************************************************************
!> \brief Deallocate all arrays pointed to by the pointers stored in the
!> real pointer array
!> \param real_arr_p ...
!> \date 2013-12-12
!> \author Lukasz Walewski
! *****************************************************************************
SUBROUTINE helium_destroy_real_arr_ptr(real_arr_p)
TYPE(real_arr_ptr), DIMENSION(:), &
POINTER :: real_arr_p
CHARACTER(len=*), PARAMETER :: routineN = 'helium_destroy_real_arr_ptr', &
routineP = moduleN//':'//routineN
INTEGER :: ip
! do not attempt deallocation on null pointer
IF (.NOT.ASSOCIATED(real_arr_p)) THEN
RETURN
END IF
! deallocate memory used by each component of the pointer array
DO ip = 1, SIZE(real_arr_p)
IF (ASSOCIATED(real_arr_p(ip)%rap)) THEN
DEALLOCATE(real_arr_p(ip)%rap)
END IF
END DO
! deallocate the memory used for pointer array itself
IF (ASSOCIATED(real_arr_p)) THEN
DEALLOCATE(real_arr_p)
END IF
RETURN
END SUBROUTINE helium_destroy_real_arr_ptr
END MODULE helium_types

198
src/motion/input_cp2k_restarts.F
View file

@ -1053,17 +1053,14 @@ CONTAINS
routineP = moduleN//':'//routineN
CHARACTER(LEN=default_string_length) :: err_str, stmp
INTEGER :: handle, i, itmp, iweight, &
INTEGER :: handle, itmp, iweight, &
msglen, nsteps, offset, reqlen
INTEGER, DIMENSION(:), POINTER :: int_msg_gather
LOGICAL :: explicit, lbf
REAL(kind=dp) :: bf, bu, rtmp
LOGICAL :: lbf
REAL(kind=dp) :: bf, bu, invproc
REAL(kind=dp), DIMENSION(3, 2) :: bg, cg, ig
REAL(kind=dp), DIMENSION(:), POINTER :: real_msg, real_msg_gather
REAL(kind=dp), DIMENSION(:, :, :), &
POINTER :: message
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: tmpsec
CALL timeset(routineN,handle)
@ -1085,15 +1082,16 @@ CONTAINS
!
! save coordinates
!
! allocate the buffer to be passed and fill it with local coords at each proc
! allocate the buffer to be passed and fill it with local coords at each
! proc
NULLIFY(real_msg)
msglen = SIZE(helium%pos)
msglen = SIZE(helium%pos(:,:,1:helium%beads))
ALLOCATE(real_msg(msglen))
real_msg(:) = PACK( helium%pos, .TRUE. )
real_msg(:) = PACK( helium%pos(:,:,1:helium%beads), .TRUE. )
! allocate the buffer for message passing
NULLIFY(real_msg_gather)
msglen = SIZE(helium%pos) * logger%para_env%num_pe
msglen = SIZE(helium%pos(:,:,1:helium%beads)) * logger%para_env%num_pe
ALLOCATE(real_msg_gather(msglen))
! pass the message from all processors to logger%para_env%source
@ -1113,7 +1111,7 @@ CONTAINS
! DEALLOCATE since this array is only used locally
DEALLOCATE(real_msg)
!
! save permutation state
!
@ -1137,6 +1135,89 @@ CONTAINS
! "The val becomes the owner of the array" - from section_vals_val_set docu
NULLIFY(int_msg_gather)
!
! save averages
!
! update the weighting factor
itmp = helium%averages_iweight
IF ( itmp .LT. 0 ) THEN
itmp = helium%current_step-helium%first_step
ELSE
itmp = itmp + helium%current_step-helium%first_step
END IF
CALL section_vals_val_set(helium%input, &
"MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", &
i_val=itmp)
! allocate the buffer for message passing
NULLIFY(real_msg_gather)
msglen = 3 * logger%para_env%num_pe
ALLOCATE(real_msg_gather(msglen))
! gather projected area from all processors
CALL mp_gather( helium%proarea%ravr, &
real_msg_gather, &
logger%para_env%source, &
logger%para_env%group )
! update it in the global input structure
CALL section_vals_val_set( helium%input, &
"MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", &
r_vals_ptr=real_msg_gather)
! allocate the buffer for message passing
NULLIFY(real_msg_gather)
msglen = 3 * logger%para_env%num_pe
ALLOCATE(real_msg_gather(msglen))
! gather projected area squared from all processors
CALL mp_gather( helium%prarea2%ravr, &
real_msg_gather, &
logger%para_env%source, &
logger%para_env%group )
! update it in the global input structure
CALL section_vals_val_set( helium%input, &
"MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", &
r_vals_ptr=real_msg_gather)
! allocate the buffer for message passing
NULLIFY(real_msg_gather)
msglen = 3 * logger%para_env%num_pe
ALLOCATE(real_msg_gather(msglen))
! gather winding number squared from all processors
CALL mp_gather( helium%wnmber2%ravr, &
real_msg_gather, &
logger%para_env%source, &
logger%para_env%group )
! update it in the global input structure
CALL section_vals_val_set( helium%input, &
"MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", &
r_vals_ptr=real_msg_gather)
! allocate the buffer for message passing
NULLIFY(real_msg_gather)
msglen = 3 * logger%para_env%num_pe
ALLOCATE(real_msg_gather(msglen))
! gather moment of inertia from all processors
CALL mp_gather( helium%mominer%ravr, &
real_msg_gather, &
logger%para_env%source, &
logger%para_env%group )
! update it in the global input structure
CALL section_vals_val_set( helium%input, &
"MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", &
r_vals_ptr=real_msg_gather)
! NULLIFY, but do not DEALLOCATE! - a new pointer to this array is silently
! assigned in section_vals_val_set - this memory will be used later on!
! "The val becomes the owner of the array" - from section_vals_val_set docu
NULLIFY(real_msg_gather)
!
! save RNG state
!
@ -1228,70 +1309,73 @@ CONTAINS
NULLIFY(real_msg_gather)
END IF
!
! save the RDFs
!
IF (helium%rdf_present) THEN
! work on the temporary array so that accumulated data remains intact
helium%rdf_inst(:,:,:) = helium%rdf_accu(:,:,:)
! average over processors / helium environments
CALL mp_sum(helium%rdf_inst,logger%para_env%group)
itmp = logger%para_env%num_pe
invproc = 1.0_dp / REAL(itmp,dp)
helium%rdf_inst(:,:,:) = helium%rdf_inst(:,:,:) * invproc
nsteps = helium%current_step-helium%first_step
helium%rdf_inst(:,:,:) = helium%rdf_inst(:,:,:) / REAL(nsteps,dp)
iweight = helium%rdf_iweight
! average over the old and the current density (observe the weights!)
helium%rdf_inst(:,:,:) = nsteps * helium%rdf_inst(:,:,:) + &
iweight * helium%rdf_rstr(:,:,:)
helium%rdf_inst(:,:,:) = helium%rdf_inst(:,:,:) / REAL(nsteps + iweight,dp)
! update in the global input structure
NULLIFY(real_msg)
msglen = SIZE(helium%rdf_inst)
ALLOCATE(real_msg(msglen))
real_msg(:) = PACK( helium%rdf_inst, .TRUE. )
CALL section_vals_val_set(&
helium%input, &
"MOTION%PINT%HELIUM%AVERAGES%RDF", &
r_vals_ptr=real_msg)
NULLIFY(real_msg)
END IF
!
! save the densities
!
IF (helium%rho_present) THEN
! make sure that the section is explicitly present in the input structure
NULLIFY(tmpsec)
tmpsec => section_vals_get_subs_vals(helium%input, &
"MOTION%PINT%HELIUM%RHO")
CALL section_vals_get(tmpsec,explicit=explicit)
IF ( .NOT. explicit ) THEN
CALL section_vals_add_values(tmpsec)
END IF
! work on the temporary array so that accumulated data remains intact
helium%rho_inst(:,:,:,:) = helium%rho_avrg(:,:,:,:)
helium%rho_inst(:,:,:,:) = helium%rho_accu(:,:,:,:)
! average all the density estimators over helium environments
DO i = 1, helium%rho_num
NULLIFY(message)
message => helium%rho_inst(i,1::1,1::1,1::1)
CALL mp_sum(message,logger%para_env%group)
END DO
! average over processors / helium environments
CALL mp_sum(helium%rho_inst,logger%para_env%group)
itmp = logger%para_env%num_pe
CPASSERT(itmp>0)
rtmp = 1.0_dp / REAL(itmp,dp)
helium%rho_inst(:,:,:,:) = helium%rho_inst(:,:,:,:) * rtmp
invproc = 1.0_dp / REAL(itmp,dp)
helium%rho_inst(:,:,:,:) = helium%rho_inst(:,:,:,:) * invproc
! average over steps performed so far in this run
nsteps = helium%current_step-helium%first_step
CPASSERT(nsteps>0)
rtmp = 1.0_dp / REAL(nsteps,dp)
helium%rho_inst(:,:,:,:) = helium%rho_inst(:,:,:,:) * rtmp
iweight = helium%rho_iweight
rtmp = 1.0_dp / REAL(nsteps+iweight,dp)
IF ( helium%rho_restart ) THEN
! average over the old and the current density (observe the weights!)
helium%rho_inst(:,:,:,:) = nsteps * helium%rho_inst(:,:,:,:) + &
iweight * helium%rho_rstr(:,:,:,:)
helium%rho_inst(:,:,:,:) = helium%rho_inst(:,:,:,:) * rtmp
END IF
helium%rho_inst(:,:,:,:) = helium%rho_inst(:,:,:,:) / REAL(nsteps,dp)
iweight = helium%averages_iweight
! average over the old and the current density (observe the weights!)
helium%rho_inst(:,:,:,:) = nsteps * helium%rho_inst(:,:,:,:) + &
iweight * helium%rho_rstr(:,:,:,:)
helium%rho_inst(:,:,:,:) = helium%rho_inst(:,:,:,:) / REAL(nsteps + iweight,dp)
! update the densities in the global input structure
NULLIFY(real_msg)
msglen = SIZE(helium%rho_inst)
ALLOCATE(real_msg(msglen))
real_msg(:) = PACK( helium%rho_inst, .TRUE. )
CALL section_vals_val_set(helium%input, &
"MOTION%PINT%HELIUM%RHO%CUBE_DATA%_DEFAULT_KEYWORD_", &
r_vals_ptr=real_msg)
CALL section_vals_val_set( &
helium%input, &
"MOTION%PINT%HELIUM%AVERAGES%RHO", &
r_vals_ptr=real_msg)
NULLIFY(real_msg)
! update the weighting factor
itmp = helium%rho_iweight
IF ( itmp .LT. 0 ) THEN
itmp = helium%current_step-helium%first_step
ELSE
itmp = itmp + helium%current_step-helium%first_step
END IF
CALL section_vals_val_set(helium%input, &
"MOTION%PINT%HELIUM%RHO%IWEIGHT", &
i_val=itmp)
END IF
CALL timestop(handle)

131
src/motion/pint_methods.F
View file

@ -26,6 +26,8 @@ MODULE pint_methods
cp_to_string
USE cp_output_handling, ONLY: cp_add_iter_level,&
cp_iterate,&
cp_print_key_finished_output,&
cp_print_key_unit_nr,&
cp_rm_iter_level
USE cp_para_types, ONLY: cp_para_env_type
USE cp_subsys_types, ONLY: cp_subsys_get,&
@ -43,6 +45,8 @@ MODULE pint_methods
gle_init,&
gle_thermo_create
USE global_types, ONLY: global_environment_type
USE helium_interactions, ONLY: helium_intpot_scan
USE helium_io, ONLY: helium_write_cubefile
USE helium_methods, ONLY: helium_create,&
helium_init,&
helium_release
@ -62,7 +66,8 @@ MODULE pint_methods
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE machine, ONLY: m_walltime
USE machine, ONLY: m_flush,&
m_walltime
USE mathconstants, ONLY: twopi
USE mathlib, ONLY: gcd
USE parallel_rng_types, ONLY: GAUSSIAN,&
@ -339,6 +344,7 @@ CONTAINS
pint_env%uv_new(pint_env%p,pint_env%ndim), &
pint_env%uf(pint_env%p,pint_env%ndim), &
pint_env%uf_h(pint_env%p,pint_env%ndim), &
pint_env%centroid(pint_env%ndim), &
pint_env%rtmp_ndim(pint_env%ndim), &
pint_env%rtmp_natom(pint_env%ndim/3), &
STAT=stat)
@ -354,6 +360,7 @@ CONTAINS
pint_env%uv_new = 0._dp
pint_env%uf = 0._dp
pint_env%uf_h = 0._dp
pint_env%centroid(:) = 0.0_dp
pint_env%rtmp_ndim= 0._dp
pint_env%rtmp_natom= 0._dp
pint_env%time_per_step = 0.0_dp
@ -610,6 +617,7 @@ CONTAINS
DEALLOCATE(pint_env%uv_new)
DEALLOCATE(pint_env%uf)
DEALLOCATE(pint_env%uf_h)
DEALLOCATE(pint_env%centroid)
DEALLOCATE(pint_env%rtmp_ndim)
DEALLOCATE(pint_env%rtmp_natom)
DEALLOCATE(pint_env%propagator)
@ -731,13 +739,13 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'do_pint_run', &
routineP = moduleN//':'//routineN
INTEGER, PARAMETER :: helium_only_mid = 1, &
solute_only_mid = 2, &
solute_with_helium_mid = 3
INTEGER, PARAMETER :: helium_only_mid = 1, int_pot_scan_mid = 4, &
solute_only_mid = 2, solute_with_helium_mid = 3
CHARACTER(len=default_string_length) :: stmp
INTEGER :: handle, mode
LOGICAL :: explicit, helium_only, &
solvent_present
int_pot_scan, solvent_present
TYPE(helium_solvent_type), POINTER :: helium
TYPE(pint_env_type), POINTER :: pint_env
TYPE(section_vals_type), POINTER :: helium_section
@ -770,13 +778,31 @@ CONTAINS
helium_only = .FALSE.
END IF
! pick the mode of operation
! check wheather to perform solute-helium interaction pot scan
IF (solvent_present) THEN
CALL section_vals_val_get(helium_section,"INTERACTION_POT_SCAN",&
l_val=int_pot_scan)
ELSE
int_pot_scan = .FALSE.
END IF
! input consistency check
IF (helium_only .AND. int_pot_scan) THEN
stmp = "Options HELIUM_ONLY and INTERACTION_POT_SCAN are exclusive"
CPABORT(stmp)
END IF
! select mode of operation
mode = 0
IF (solvent_present) THEN
IF (helium_only) THEN
mode = helium_only_mid
ELSE
mode = solute_with_helium_mid
IF (int_pot_scan) THEN
mode = int_pot_scan_mid
ELSE
mode = solute_with_helium_mid
END IF
END IF
ELSE
mode = solute_only_mid
@ -799,6 +825,17 @@ CONTAINS
CALL pint_do_run(pint_env,globenv)
CALL pint_release(pint_env)
CASE (int_pot_scan_mid)
CALL pint_create(pint_env,input,input_declaration,para_env)
! TODO only initialization of positions is necessary, but rep_env_calc_e_f called
! from within pint_init_f does something to the replica environments which can not be
! avoided (has something to do with f_env_add_defaults) so leaving for now..
CALL pint_init(pint_env)
CALL helium_create(helium,input,solute=pint_env)
CALL pint_run_scan(pint_env, globenv, helium)
CALL helium_release(helium)
CALL pint_release(pint_env)
CASE (solute_with_helium_mid)
CALL pint_create(pint_env,input,input_declaration,para_env)
! init pint wihtout helium forces (they are not yet initialized)
@ -1491,6 +1528,19 @@ CONTAINS
CALL cp_iterate(pint_env%logger%iter_info,iter_nr=pint_env%first_step)
pint_env%iter = pint_env%first_step
IF (PRESENT(helium)) THEN
IF (ASSOCIATED(helium)) THEN
! set the properties accumulated over the whole MC process to 0
helium%proarea%accu(:) = 0.0_dp
helium%prarea2%accu(:) = 0.0_dp
helium%wnmber2%accu(:) = 0.0_dp
helium%mominer%accu(:) = 0.0_dp
IF (helium%rho_present) THEN
helium%rho_accu(:,:,:,:) = 0.0_dp
END IF
END IF
END IF
! write the properties at 0-th step
CALL pint_calc_energy(pint_env)
CALL pint_write_ener(pint_env)
@ -1539,6 +1589,73 @@ CONTAINS
RETURN
END SUBROUTINE pint_do_run
! ***************************************************************************
!> \brief Performs a scan of the helium-solute interaction energy
!> \author Lukasz Walewski
!> \date 2013-11-26
! ***************************************************************************
! *****************************************************************************
!> \brief ...
!> \param pint_env ...
!> \param globenv ...
!> \param helium ...
! *****************************************************************************
SUBROUTINE pint_run_scan(pint_env, globenv, helium)
TYPE(pint_env_type), POINTER :: pint_env
TYPE(global_environment_type), POINTER :: globenv
TYPE(helium_solvent_type), POINTER :: helium
CHARACTER(len=*), PARAMETER :: routineN = 'pint_run_scan', &
routineP = moduleN//':'//routineN
CHARACTER(len=default_string_length) :: comment
INTEGER :: unit_nr
REAL(kind=dp), DIMENSION(:, :, :), &
POINTER :: DATA
TYPE(section_vals_type), POINTER :: print_key
NULLIFY(pint_env%logger, print_key)
pint_env%logger => cp_get_default_logger()
print_key => section_vals_get_subs_vals(helium%input, &
"MOTION%PINT%HELIUM%PRINT%RHO")
! perform the actual scan wrt the COM of the solute
helium%rho_inst(1,:,:,:) = 0.0_dp
CALL helium_intpot_scan(pint_env, helium)
! output the interaction potential into a cubefile
IF (pint_env%logger%para_env%ionode) THEN
unit_nr=cp_print_key_unit_nr( &
pint_env%logger, &
print_key, &
middle_name="helium-pot", &
extension=".cube", &
file_position="REWIND", &
do_backup=.FALSE.)
comment = "Solute - helium interaction potential"
NULLIFY(DATA)
DATA => helium%rho_inst(1,:,:,:)
CALL helium_write_cubefile( &
unit_nr, &
comment, &
helium%center-0.5_dp*(helium%rho_maxr-helium%rho_delr), &
helium%rho_delr, &
helium%rho_nbin, &
DATA )
CALL m_flush(unit_nr)
CALL cp_print_key_finished_output(unit_nr,pint_env%logger,print_key)
END IF
! output solute positions
CALL pint_write_centroids(pint_env)
CALL pint_write_trajectory(pint_env)
RETURN
END SUBROUTINE pint_run_scan
! ***************************************************************************
!> \brief Does an MD step (and nrespa harmonic evaluations)

32
src/motion/pint_public.F
View file

@ -27,6 +27,7 @@ MODULE pint_public
PUBLIC :: pint_com_pos
PUBLIC :: pint_levy_walk
PUBLIC :: pint_calc_centroid
CONTAINS
@ -268,4 +269,35 @@ MODULE pint_public
RETURN
END SUBROUTINE pint_levy_walk
! ***************************************************************************
!> \brief Calculate the centroid
!> \param pint_env path integral environment
!> \date 2013-02-10
!> \author lwalewski
! ***************************************************************************
! *****************************************************************************
!> \brief ...
!> \param pint_env ...
! *****************************************************************************
SUBROUTINE pint_calc_centroid(pint_env)
TYPE(pint_env_type), POINTER :: pint_env
CHARACTER(len=*), PARAMETER :: routineN = 'pint_calc_centroid', &
routineP = moduleN//':'//routineN
INTEGER :: ia, ib
REAL(KIND=dp) :: invp
invp = 1.0_dp / pint_env%p
pint_env%centroid(:) = 0.0_dp
DO ia = 1, pint_env%ndim
DO ib = 1, pint_env%p
pint_env%centroid(ia) = pint_env%centroid(ia) + pint_env%x(ib,ia)
END DO
pint_env%centroid(ia) = pint_env%centroid(ia) * invp
END DO
END SUBROUTINE pint_calc_centroid
END MODULE pint_public

2
src/motion/pint_types.F
View file

@ -122,7 +122,7 @@ MODULE pint_types
REAL(KIND=dp), DIMENSION(:), POINTER :: mass,e_pot_bead
REAL(KIND=dp), DIMENSION(:,:), POINTER :: x,v,f,mass_beads,&
mass_fict,ux,ux_t,uv,uv_t,uv_new,uf,uf_h, external_f
REAL(KIND=dp), DIMENSION(:), POINTER :: centroid
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: tx,tv,tv_t,tv_old,tv_new,tf
REAL(KIND=dp), DIMENSION(:), POINTER :: Q ! dim p, make it (p,ndim)?
REAL(KIND=dp), DIMENSION(:), POINTER :: rtmp_ndim, rtmp_natom

5
src/start/input_cp2k.F
View file

@ -599,6 +599,11 @@ CONTAINS
default_l_val=.FALSE.)
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="RESTART_HELIUM_AVERAGES",&
description="Restarts average properties from PINT%HELIUM%AVERAGES.",&
type_of_var=logical_t,lone_keyword_l_val=.TRUE.,default_l_val=.FALSE.)
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_ext_restart_section

372
src/start/input_cp2k_motion.F
View file

@ -27,9 +27,14 @@ MODULE input_cp2k_motion
default_minimization_method_id, default_ts_method_id, do_mc_gemc_npt, &
do_mc_gemc_nvt, do_mc_traditional, do_mc_virial, fmt_id_pdb, &
fmt_id_xyz, gaussian, helium_cell_shape_cube, &
helium_cell_shape_octahedron, integrate_exact, integrate_numeric, &
ls_2pnt, ls_3pnt, ls_fit, ls_gold, ls_none, matrix_init_cholesky, &
matrix_init_diagonal, numerical, propagator_pimd, propagator_rpmd, &
helium_cell_shape_octahedron, helium_forces_average, &
helium_forces_last, helium_mdist_exponential, helium_mdist_gaussian, &
helium_mdist_linear, helium_mdist_quadratic, helium_mdist_singlev, &
helium_mdist_uniform, helium_sampling_ceperley, &
helium_solute_intpot_mwater, helium_solute_intpot_none, &
integrate_exact, integrate_numeric, ls_2pnt, ls_3pnt, ls_fit, ls_gold, &
ls_none, matrix_init_cholesky, matrix_init_diagonal, numerical, &
perm_cycle, perm_plain, propagator_pimd, propagator_rpmd, &
transformation_normal, transformation_stage
USE input_cp2k_constraints, ONLY: create_constraint_section
USE input_cp2k_free_energy, ONLY: create_fe_section
@ -51,6 +56,7 @@ MODULE input_cp2k_motion
section_release,&
section_type
USE input_val_types, ONLY: integer_t,&
logical_t,&
real_t
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
@ -1922,15 +1928,14 @@ CONTAINS
routineP = moduleN//':'//routineN
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: print_key, subsection, &
subsubsection
TYPE(section_type), POINTER :: print_key, subsection
CPASSERT(.NOT.ASSOCIATED(section))
CALL section_create(section,name="HELIUM",&
description="The section that controls optional helium solvent"// &
" environment (highly experimental, not for general use yet)",&
n_keywords=11, n_subsections=4, repeats=.FALSE.)
n_keywords=28, n_subsections=9, repeats=.FALSE.)
NULLIFY(keyword)
CALL keyword_create(keyword, name="_SECTION_PARAMETERS_",&
@ -1947,6 +1952,14 @@ CONTAINS
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="INTERACTION_POT_SCAN",&
description="Scan solute-helium interaction potential, "//&
"cubefile parameters set in subsection RHO",&
repeats=.FALSE., default_l_val=.FALSE.,&
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="NUM_ENV", &
description="Number of independent helium environments"// &
" (only for restarts, do not set explicitly)",&
@ -1960,6 +1973,29 @@ CONTAINS
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="GET_FORCES",&
description="Get average MC forces or last MC forces to propagate MD",&
usage="GET_FORCES (AVERAGE|LAST)",&
default_i_val=helium_forces_average,&
enum_c_vals=s2a("AVERAGE","LAST"),&
enum_i_vals=(/helium_forces_average,helium_forces_last/))
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="SOLUTE_INTERACTION", &
description="Interaction potential between helium and the solute",&
usage="SOLUTE_INTERACTION (NONE | MWATER)",&
default_i_val=helium_solute_intpot_none,&
enum_c_vals=s2a("NONE", "MWATER"),&
enum_i_vals=(/&
helium_solute_intpot_none, &
helium_solute_intpot_mwater/), &
enum_desc=s2a(&
"No interaction with solute", &
"Test interaction with wrong Water"))
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="NATOMS", &
description="Number of helium atoms",&
repeats=.FALSE., default_i_val=64)
@ -1972,15 +2008,19 @@ CONTAINS
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="INOROT", &
description="Number of MC iterations at the same time slice(s)",&
repeats=.FALSE., default_i_val=10000)
CALL keyword_create(keyword, name="N_INNER", &
variants=s2a("INOROT"),&
description="Number of MC iterations at the same time slice(s) "//&
"(number of inner MC loop iterations)",&
repeats=.FALSE., default_i_val=6600)
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="IROT", &
description="how often to reselect the time slice(s) to work on",&
repeats=.FALSE., default_i_val=10000)
CALL keyword_create(keyword, name="N_OUTER", &
variants=s2a("IROT"),&
description="how often to reselect the time slice(s) to work on "//&
"(number of outer MC loop iterations)",&
repeats=.FALSE., default_i_val=300)
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
@ -1997,12 +2037,42 @@ CONTAINS
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="SAMPLING_METHOD",&
description="Choose between Ceperley or other not implemented algorithms",&
usage="SAMPLING_METHOD (CEPERLEY)",&
default_i_val=helium_sampling_ceperley,&
enum_c_vals=s2a("CEPERLEY"),&
enum_i_vals=(/helium_sampling_ceperley/))
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
NULLIFY(subsection)
CALL section_create(subsection,name="M-SAMPLING",&
description="Permutation cycle length sampling settings",&
n_keywords=2, n_subsections=0, repeats=.FALSE.)
n_keywords=3, n_subsections=0, repeats=.FALSE.)
CALL keyword_create(keyword, name="DISTRIBUTION-TYPE", &
description="Distribution from which the cycle length m is sampled",&
usage="DISTRIBUTION-TYPE (SINGLEV|UNIFORM|LINEAR|QUADRATIC|EXPONENTIAL|GAUSSIAN)",&
default_i_val=helium_mdist_uniform,&
enum_c_vals=s2a(&
"SINGLEV",&
"UNIFORM",&
"LINEAR",&
"QUADRATIC",&
"EXPONENTIAL",&
"GAUSSIAN"),&
enum_i_vals=(/&
helium_mdist_singlev,&
helium_mdist_uniform,&
helium_mdist_linear,&
helium_mdist_quadratic,&
helium_mdist_exponential,&
helium_mdist_gaussian/))
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="M-VALUE", &
description="Value of m treated in a special way",&
description="Value of m treated in a special way "//&
"(specific behavior depends on the distribution type chosen)",&
repeats=.FALSE.,&
default_i_val=1)
CALL section_add_keyword(subsection,keyword)
@ -2039,6 +2109,14 @@ CONTAINS
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="DROPLET_RADIUS", &
description="Reject a move if any of the new positions does not lie within"//&
" this range from the center of gravity",&
repeats=.FALSE.,type_of_var=real_t,default_r_val=HUGE(1.0_dp),&
unit_str="angstrom")
CALL section_add_keyword(section,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="DENSITY", &
description="trial density of helium for determining the helium "//&
"box size",&
@ -2062,27 +2140,42 @@ CONTAINS
CALL keyword_release(keyword)
CALL section_create(subsection,name="RDF",&
description="Radial distribution function generation settings",&
n_keywords=2, n_subsections=0, repeats=.FALSE.)
description="Radial distribution settings",&
n_keywords=4, n_subsections=0, repeats=.FALSE.)
CALL keyword_create(keyword, name="_SECTION_PARAMETERS_",&
description="Whether or not to actually calculate this property",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="MAXR", &
description="Maximum RDF range, defaults to unit cell size",&
repeats=.FALSE.,type_of_var=real_t,&
unit_str="angstrom")
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="NBIN", &
description="Number of bins",&
repeats=.FALSE.,&
default_i_val=700)
default_i_val=250)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="CENTERS_FILE_NAME",&
description="Calculate RDFs wrt centers read from a file",&
repeats=.FALSE., default_lc_val="centers.xyz")
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section,subsection)
CALL section_release(subsection)
NULLIFY(subsection)
CALL section_create(subsection,name="RHO",&
description="Density distribution settings",&
n_keywords=2, n_subsections=0, repeats=.FALSE.)
description="Spatial distribution settings",&
n_keywords=10, n_subsections=0, repeats=.FALSE.)
CALL keyword_create(keyword, name="_SECTION_PARAMETERS_",&
description="Whether or not to actually calculate densities "//&
"(requires significant amount of memory, depending on the value of NBIN)",&
@ -2090,33 +2183,80 @@ CONTAINS
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="NBIN", &
description="Number of bins",&
description="Number of grid points in each direction for density binning",&
repeats=.FALSE.,&
default_i_val=100)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="IWEIGHT", &
description="Weight the restarted density should be given " // &
"(number of MC steps used to average the restarted density, " // &
"negative value - the same weight as the run-time density, " // &
"usually should not be changed)",&
repeats=.FALSE.,&
default_i_val=-1)
!
CALL keyword_create(keyword, name="MIN_CYCLE_LENGTHS_WDG",&
description="Density of winding paths "// &
"not shorter than the given length", &
repeats=.FALSE.,usage="<INT> <INT> .. <INT>",&
type_of_var=integer_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
NULLIFY(subsubsection)
CALL section_create(subsubsection,name="CUBE_DATA",&
description="Density data used for restarts",&
n_keywords=1, n_subsections=0, repeats=.FALSE.)
CALL keyword_create(keyword, name="_DEFAULT_KEYWORD_",&
description="Cubefile data", &
repeats=.TRUE., usage="{Real} ...",&
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(subsubsection,keyword)
!
CALL keyword_create(keyword, name="MIN_CYCLE_LENGTHS_NON",&
description="Density of non-winding paths "// &
"not shorter than the given length", &
repeats=.FALSE.,usage="<INT> <INT> .. <INT>",&
type_of_var=integer_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="MIN_CYCLE_LENGTHS_ALL",&
description="Density of all paths "// &
"not shorter than the given length", &
repeats=.FALSE.,usage="<INT> <INT> .. <INT>",&
type_of_var=integer_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="ATOM_NUMBER",&
description="Atom number density",&
repeats=.FALSE.,&
type_of_var=logical_t,&
default_l_val=.TRUE.,&
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="PROJECTED_AREA_2",&
description="Projected area squared density, A*A(r)",&
repeats=.FALSE.,&
type_of_var=logical_t,&
default_l_val=.FALSE.,&
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="WINDING_NUMBER_2",&
description="Winding number squared density, W*W(r)",&
repeats=.FALSE.,&
type_of_var=logical_t,&
default_l_val=.FALSE.,&
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="WINDING_CYCLE_2",&
description="Winding number squared density, W^2(r)",&
repeats=.FALSE.,&
type_of_var=logical_t,&
default_l_val=.FALSE.,&
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
!
CALL keyword_create(keyword, name="MOMENT_OF_INERTIA",&
description="Moment of inertia density",&
repeats=.FALSE.,&
type_of_var=logical_t,&
default_l_val=.FALSE.,&
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(subsection,subsubsection)
CALL section_release(subsubsection)
CALL section_add_subsection(section,subsection)
CALL section_release(subsection)
@ -2131,21 +2271,70 @@ CONTAINS
n_keywords=1, n_subsections=0, repeats=.FALSE.)
CALL keyword_create(keyword, name="_DEFAULT_KEYWORD_",&
description="Specify particle index permutation for every "// &
"helium atom",repeats=.TRUE.,usage="{Integer} ...",&
"helium atom",repeats=.TRUE.,usage="<INT> <INT> .. <INT>",&
type_of_var=integer_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section,subsection)
CALL section_release(subsection)
CALL section_create(subsection,name="AVERAGES",&
description="Average properties (used for restarts)",&
n_keywords=7, n_subsections=0, repeats=.FALSE.)
CALL keyword_create(keyword, name="PROJECTED_AREA",&
description="Projected area vector for all environments", &
repeats=.TRUE.,usage="<REAL> <REAL> .. <REAL>",&
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="PROJECTED_AREA_2",&
description="Projected area vector squared for all environments", &
repeats=.TRUE.,usage="<REAL> <REAL> .. <REAL>",&
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="WINDING_NUMBER_2",&
description="Winding number vector squared for all environments", &
repeats=.TRUE.,usage="<REAL> <REAL> .. <REAL>",&
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="MOMENT_OF_INERTIA",&
description="Moment of inertia vector for all environments", &
repeats=.TRUE.,usage="<REAL> <REAL> .. <REAL>",&
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="RDF",&
description="Radial distributions averaged over all environments", &
repeats=.TRUE.,usage="<REAL> <REAL> .. <REAL>",&
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="RHO",&
description="Spatial distributions averaged over all environments", &
repeats=.TRUE.,usage="<REAL> <REAL> .. <REAL>",&
type_of_var=real_t, n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="IWEIGHT", &
description="Weight for the restarted quantities " // &
"(number of MC steps used to calculate the accumulated averages)", &
repeats=.FALSE.,&
default_i_val=0)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section,subsection)
CALL section_release(subsection)
CALL section_create(subsection,name="FORCE",&
description="Forces exerted by the helium on the solute system"//&
" (used for restarts)",&
n_keywords=0, n_subsections=0, repeats=.FALSE.)
n_keywords=1, n_subsections=0, repeats=.FALSE.)
CALL keyword_create(keyword, name="_DEFAULT_KEYWORD_", &
description="Number of real values should be 3 * "//&
"<num_solute_atoms> * <num_solute_beads>", repeats=.TRUE., &
usage="{Real} ...", type_of_var=real_t, &
usage="<REAL> <REAL> .. <REAL>", type_of_var=real_t, &
n_var=-1)
CALL section_add_keyword(subsection,keyword)
CALL keyword_release(keyword)
@ -2167,7 +2356,7 @@ CONTAINS
CALL section_create(subsection,name="PRINT",&
description="The section that controls the output of the helium code",&
n_keywords=1, n_subsections=0, repeats=.FALSE.)
n_keywords=15, n_subsections=0, repeats=.FALSE.)
! *************************************************************************
!> Printkeys for properites output
@ -2180,38 +2369,35 @@ CONTAINS
! Properties printed at LOW print level
!
CALL cp_print_key_section_create(print_key,"ENERGY",&
description="Controls the output of the helium energies"//&
description="Controls the output of helium energies"//&
" (averaged over MC step)", &
print_level=low_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key,"SDENSITY",&
description="Controls the output of the helium superfluid density",&
!
CALL cp_print_key_section_create(print_key,"PROJECTED_AREA_2_AVG",&
description="Controls the output of the average projected area squared vector",&
print_level=low_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
!
CALL cp_print_key_section_create(print_key,"WINDING_NUMBER_2_AVG",&
description="Controls the output of the average winding number vector squared",&
print_level=low_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
!
CALL cp_print_key_section_create(print_key,"MOMENT_OF_INERTIA_AVG",&
description="Controls the output of the average moment of inertia vector",&
print_level=low_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
! Properties printed at MEDIUM print level
!
CALL cp_print_key_section_create(print_key,"COORDINATES",&
description="Controls the output of helium coordinates",&
print_level=medium_print_level,common_iter_levels=1)
CALL keyword_create(keyword, name="FORMAT",&
description="Output file format for the coordinates",&
usage="FORMAT (PDB|XYZ)",&
default_i_val=fmt_id_pdb,&
enum_c_vals=s2a("PDB","XYZ"),&
enum_i_vals= (/fmt_id_pdb,fmt_id_xyz/),&
enum_desc=s2a( "Bead coordinates and connectivity is written in PDB format",&
"Only bead coordinates are written in XYZ format"))
CALL section_add_keyword(print_key,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key,"RDF",&
description="Controls the output of the helium radial distribution function",&
description="Controls the output of helium radial distribution functions",&
print_level=medium_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
@ -2230,7 +2416,25 @@ CONTAINS
CALL keyword_release(keyword)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
!
CALL cp_print_key_section_create(print_key,"PROJECTED_AREA",&
description="Controls the output of the projected area vector",&
print_level=medium_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
!
CALL cp_print_key_section_create(print_key,"WINDING_NUMBER",&
description="Controls the output of the winding number vector",&
print_level=medium_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
!
CALL cp_print_key_section_create(print_key,"MOMENT_OF_INERTIA",&
description="Controls the output of the moment of inertia vector",&
print_level=medium_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
!
CALL cp_print_key_section_create(print_key,"PLENGTH",&
description="Controls the output of the helium permutation length",&
print_level=medium_print_level,common_iter_levels=1)
@ -2239,15 +2443,37 @@ CONTAINS
! Properties printed at HIGH print level
!
CALL cp_print_key_section_create(print_key,"ACCEPTS",&
description="Controls the output of the helium acceptance data",&
CALL cp_print_key_section_create(print_key,"COORDINATES",&
description="Controls the output of helium coordinates",&
print_level=high_print_level,common_iter_levels=1)
CALL keyword_create(keyword, name="FORMAT",&
description="Output file format for the coordinates",&
usage="FORMAT (PDB|XYZ)",&
default_i_val=fmt_id_pdb,&
enum_c_vals=s2a("PDB","XYZ"),&
enum_i_vals= (/fmt_id_pdb,fmt_id_xyz/),&
enum_desc=s2a( "Bead coordinates and connectivity is written in PDB format",&
"Only bead coordinates are written in XYZ format"))
CALL section_add_keyword(print_key,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
!
CALL cp_print_key_section_create(print_key,"PERM",&
description="Controls the output of the helium permutation state",&
print_level=high_print_level,common_iter_levels=1)
CALL keyword_create(keyword, name="FORMAT",&
description="Output format for the permutation",&
usage="FORMAT (CYCLE|PLAIN)",&
default_i_val=perm_cycle,&
enum_c_vals=s2a("CYCLE","PLAIN"),&
enum_i_vals= (/perm_cycle,perm_plain/),&
enum_desc=s2a( &
"Cycle notation with winding cycles enclosed"//&
" in '[...]' and non-winding ones enclosed in '(...)'",&
"Plain permutation output, i.e. P(1) ... P(N)"))
CALL section_add_keyword(print_key,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
@ -2259,14 +2485,14 @@ CONTAINS
! Properties printed at DEBUG print level
!
CALL cp_print_key_section_create(print_key,"FORCES_INST",&
description="Controls the output of the instantaneous helium forces on the solute",&
CALL cp_print_key_section_create(print_key,"ACCEPTS",&
description="Controls the output of the helium acceptance data",&
print_level=debug_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key,"WNUMBER",&
description="Controls the output of the helium winding number",&
!
CALL cp_print_key_section_create(print_key,"FORCES_INST",&
description="Controls the output of the instantaneous helium forces on the solute",&
print_level=debug_print_level,common_iter_levels=1)
CALL section_add_subsection(subsection,print_key)
CALL section_release(print_key)

1
tests/Pimd/TEST_FILES
View file

@ -12,6 +12,7 @@ he32_only.inp 40
he32_only_restart.inp 40 2e-14
water-in-helium.inp 9 7e-14
water-in-helium-restart.inp 9 7e-14
water-in-helium-lastforce.inp 9 7e-14
w512_pint_nose.inp 9 8e-14 9.5780711530210407
w512_pint_gle.inp 9 3e-13 -3.5646497002020765
w512_pint_pile.inp 9 3e-12 0.70661010702892924

67
tests/Pimd/water-in-helium-lastforce.inp Normal file
View file

@ -0,0 +1,67 @@
&GLOBAL
PROJECT_NAME water-in-helium-lastforce
RUN_TYPE PINT
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&PINT
P 4
PROC_PER_REPLICA 1
NUM_STEPS 5
DT 0.5
NRESPA 2
TEMP 7.5
T_TOL 0.0
TRANSFORMATION NORMAL
&NOSE
NNOS 3
&END NOSE
&HELIUM
SOLUTE_INTERACTION MWATER
POTENTIAL_FILE_NAME lj-test-pot.dat
GET_FORCES LAST
NATOMS 32
NBEADS 16
INOROT 200
IROT 100
BISECTION 8
MAX_PERM_CYCLE 4
PERIODIC T
CELL_SHAPE OCTAHEDRON
&END HELIUM
&END PINT
&END MOTION
&FORCE_EVAL
METHOD FIST
&MM
&POISSON
&EWALD
GMAX 25
&END EWALD
&END POISSON
&FORCEFIELD
PARM_FILE_NAME ../Fist/sample_pot/water.pot
PARMTYPE CHM
&CHARGE
ATOM OT
CHARGE -0.8476
&END CHARGE
&CHARGE
ATOM HT
CHARGE 0.4238
&END CHARGE
&END FORCEFIELD
&END MM
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME ../Fist/sample_pdb/water_1.pdb
COORD_FILE_FORMAT PDB
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

1
tests/Pimd/water-in-helium-restart.inp
View file

@ -20,6 +20,7 @@
NNOS 3
&END NOSE
&HELIUM
SOLUTE_INTERACTION MWATER
POTENTIAL_FILE_NAME lj-test-pot.dat
NATOMS 32
NBEADS 16

1
tests/Pimd/water-in-helium.inp
View file

@ -18,6 +18,7 @@
NNOS 3
&END NOSE
&HELIUM
SOLUTE_INTERACTION MWATER
POTENTIAL_FILE_NAME lj-test-pot.dat
NATOMS 32
NBEADS 16