1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Methods dealing with helium_solvent_type 8!> \author Lukasz Walewski 9!> \date 2009-06-10 10! ************************************************************************************************** 11MODULE helium_methods 12 13 USE atomic_kind_types, ONLY: get_atomic_kind 14 USE bibliography, ONLY: Walewski2014,& 15 cite_reference 16 USE cell_types, ONLY: cell_type,& 17 get_cell 18 USE cp_files, ONLY: close_file,& 19 open_file 20 USE cp_log_handling, ONLY: cp_get_default_logger,& 21 cp_logger_type 22 USE cp_output_handling, ONLY: cp_printkey_is_on 23 USE cp_subsys_types, ONLY: cp_subsys_get,& 24 cp_subsys_type 25 USE f77_interface, ONLY: f_env_add_defaults,& 26 f_env_rm_defaults,& 27 f_env_type 28 USE force_env_types, ONLY: force_env_get 29 USE helium_common, ONLY: helium_com,& 30 helium_pbc 31 USE helium_interactions, ONLY: helium_vij 32 USE helium_io, ONLY: helium_write_line,& 33 helium_write_setup 34 USE helium_sampling, ONLY: helium_sample 35 USE helium_types, ONLY: helium_solvent_p_type,& 36 helium_solvent_type,& 37 rho_atom_number,& 38 rho_moment_of_inertia,& 39 rho_num,& 40 rho_projected_area,& 41 rho_winding_cycle,& 42 rho_winding_number 43 USE input_constants, ONLY: helium_cell_shape_cube,& 44 helium_cell_shape_octahedron,& 45 helium_sampling_ceperley,& 46 helium_sampling_worm,& 47 helium_solute_intpot_none 48 USE input_section_types, ONLY: section_vals_get,& 49 section_vals_get_subs_vals,& 50 section_vals_type,& 51 section_vals_val_get,& 52 section_vals_val_set 53 USE kinds, ONLY: default_path_length,& 54 default_string_length,& 55 dp,& 56 max_line_length 57 USE mathconstants, ONLY: pi,& 58 twopi 59 USE message_passing, ONLY: mp_allgather,& 60 mp_bcast,& 61 mp_comm_free,& 62 mp_comm_split_direct 63 USE parallel_rng_types, ONLY: GAUSSIAN,& 64 UNIFORM,& 65 create_rng_stream,& 66 delete_rng_stream,& 67 next_random_number,& 68 rng_stream_p_type,& 69 rng_stream_type,& 70 set_rng_stream 71 USE particle_list_types, ONLY: particle_list_type 72 USE physcon, ONLY: a_mass,& 73 angstrom,& 74 boltzmann,& 75 h_bar,& 76 kelvin,& 77 massunit 78 USE pint_public, ONLY: pint_com_pos 79 USE pint_types, ONLY: pint_env_type 80 USE splines_methods, ONLY: init_spline,& 81 init_splinexy 82 USE splines_types, ONLY: spline_data_create,& 83 spline_data_release 84#include "../base/base_uses.f90" 85 86 IMPLICIT NONE 87 88 PRIVATE 89 90 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE. 91 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_methods' 92 INTEGER, SAVE, PRIVATE :: last_helium_id = 0 93 94 PUBLIC :: helium_create 95 PUBLIC :: helium_init 96 PUBLIC :: helium_release 97 98CONTAINS 99 100! *************************************************************************** 101!> \brief Data-structure that holds all needed information about 102!> (superfluid) helium solvent 103!> \param helium_env ... 104!> \param input ... 105!> \param solute ... 106!> \par History 107!> 2016-07-14 Modified to work with independent helium_env [cschran] 108!> \author hforbert 109! ************************************************************************************************** 110 SUBROUTINE helium_create(helium_env, input, solute) 111 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 112 TYPE(section_vals_type), POINTER :: input 113 TYPE(pint_env_type), OPTIONAL, POINTER :: solute 114 115 CHARACTER(len=*), PARAMETER :: routineN = 'helium_create', routineP = moduleN//':'//routineN 116 117 CHARACTER(len=default_path_length) :: msg_str, potential_file_name 118 INTEGER :: color_sub, handle, i, input_unit, isize, & 119 itmp, j, k, mepos, new_comm, nlines, & 120 ntab, num_env, pdx 121 INTEGER, DIMENSION(:), POINTER :: env_all 122 LOGICAL :: expl_cell, expl_dens, expl_nats, & 123 expl_pot, explicit, ltmp 124 REAL(KIND=dp) :: cgeof, dx, he_mass, mHe, rtmp, T, tau, & 125 tcheck, x1, x_spline 126 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pot_transfer 127 TYPE(cp_logger_type), POINTER :: logger 128 TYPE(section_vals_type), POINTER :: helium_section, input_worm 129 130 CALL timeset(routineN, handle) 131 132 CALL cite_reference(Walewski2014) 133 NULLIFY (logger) 134 logger => cp_get_default_logger() 135 136 NULLIFY (helium_section) 137 helium_section => section_vals_get_subs_vals(input, & 138 "MOTION%PINT%HELIUM") 139 CALL section_vals_get(helium_section, explicit=explicit) 140 CPASSERT(explicit) 141 142 ! get number of environments 143 CALL section_vals_val_get(helium_section, "NUM_ENV", & 144 explicit=explicit) 145 IF (explicit) THEN 146 CALL section_vals_val_get(helium_section, "NUM_ENV", & 147 i_val=num_env) 148 ELSE 149 num_env = logger%para_env%num_pe 150 END IF 151 CPASSERT(num_env >= 0) 152 IF (num_env .NE. logger%para_env%num_pe) THEN 153 msg_str = "NUM_ENV not equal to number of processors" 154 CPWARN(msg_str) 155 END IF 156 ! calculate number of tasks for each processor 157 mepos = num_env/logger%para_env%num_pe & 158 + MIN(MOD(num_env, logger%para_env%num_pe)/(logger%para_env%mepos + 1), 1) 159 ! gather result 160 NULLIFY (env_all) 161 ALLOCATE (env_all(logger%para_env%num_pe)) 162 env_all(:) = 0 163 CALL mp_allgather(mepos, env_all, logger%para_env%group) 164 165 ! create new communicator for processors with helium_env 166 IF (mepos == 0) THEN 167 color_sub = 0 168 ELSE 169 color_sub = 1 170 END IF 171 CALL mp_comm_split_direct(logger%para_env%group, new_comm, color_sub) 172 ! release new_comm for processors without helium_env 173 IF (mepos == 0) THEN 174 CALL mp_comm_free(new_comm) 175 END IF 176 177 NULLIFY (helium_env) 178 IF (mepos .GT. 0) THEN 179 ALLOCATE (helium_env(mepos)) 180 DO k = 1, mepos 181 !NULLIFY (helium_env%para_env) 182 !CALL cp_para_env_create(helium_env%para_env, new_comm) 183 helium_env(k)%comm = new_comm 184 NULLIFY (helium_env(k)%env_all) 185 helium_env(k)%env_all => env_all 186 ALLOCATE (helium_env(k)%helium) 187 NULLIFY (helium_env(k)%helium%input) 188 helium_env(k)%helium%input => input 189 helium_env(k)%helium%num_env = num_env 190 END DO 191 ! RNG state create & init 192 CALL helium_rng_init(helium_env) 193 194 DO k = 1, mepos 195 NULLIFY (helium_env(k)%helium%ptable, helium_env(k)%helium%permutation, & 196 helium_env(k)%helium%iperm, & 197 helium_env(k)%helium%itmp_atoms_1d, & 198 helium_env(k)%helium%ltmp_atoms_1d, & 199 helium_env(k)%helium%itmp_atoms_np_1d, & 200 helium_env(k)%helium%pos, helium_env(k)%helium%work, & 201 helium_env(k)%helium%force_avrg, & 202 helium_env(k)%helium%force_inst, & 203 helium_env(k)%helium%rtmp_3_np_1d, & 204 helium_env(k)%helium%rtmp_p_ndim_1d, & 205 helium_env(k)%helium%rtmp_p_ndim_np_1d, & 206 helium_env(k)%helium%rtmp_3_atoms_beads_1d, & 207 helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, & 208 helium_env(k)%helium%rtmp_p_ndim_2d, & 209 helium_env(k)%helium%ltmp_3_atoms_beads_3d, & 210 helium_env(k)%helium%tmatrix, helium_env(k)%helium%pmatrix, & 211 helium_env(k)%helium%nmatrix, helium_env(k)%helium%ipmatrix, & 212 helium_env(k)%helium%u0, helium_env(k)%helium%e0, & 213 helium_env(k)%helium%uoffdiag, helium_env(k)%helium%eoffdiag, & 214 helium_env(k)%helium%vij, & 215 helium_env(k)%helium%rdf_inst, & 216 helium_env(k)%helium%plength_avrg, & 217 helium_env(k)%helium%plength_inst, & 218 helium_env(k)%helium%atom_plength, & 219 helium_env(k)%helium%ename & 220 ) 221 222 helium_env(k)%helium%ref_count = 1 223 last_helium_id = last_helium_id + 1 224 helium_env(k)%helium%id_nr = last_helium_id 225 helium_env(k)%helium%accepts = 0 226 helium_env(k)%helium%relrot = 0 227 228 ! check if solute is present in our simulation 229 helium_env(k)%helium%solute_present = .FALSE. 230 helium_env(k)%helium%solute_atoms = 0 231 helium_env(k)%helium%solute_beads = 0 232 CALL section_vals_val_get( & 233 helium_section, & 234 "HELIUM_ONLY", & 235 l_val=ltmp) 236 IF (.NOT. ltmp) THEN 237 IF (PRESENT(solute)) THEN 238 IF (ASSOCIATED(solute)) THEN 239 helium_env(k)%helium%solute_present = .TRUE. 240 helium_env(k)%helium%solute_atoms = solute%ndim/3 241 helium_env(k)%helium%solute_beads = solute%p 242 END IF 243 END IF 244 END IF 245 246 CALL section_vals_val_get(helium_section, "NBEADS", & 247 i_val=helium_env(k)%helium%beads) 248 CALL section_vals_val_get(helium_section, "INOROT", & 249 i_val=helium_env(k)%helium%iter_norot) 250 CALL section_vals_val_get(helium_section, "IROT", & 251 i_val=helium_env(k)%helium%iter_rot) 252 253 ! get number of steps and current step number from PINT 254 CALL section_vals_val_get(input, "MOTION%PINT%ITERATION", & 255 i_val=itmp) 256 helium_env(k)%helium%first_step = itmp 257 CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", & 258 explicit=explicit) 259 IF (explicit) THEN 260 CALL section_vals_val_get(input, "MOTION%PINT%MAX_STEP", & 261 i_val=itmp) 262 helium_env(k)%helium%last_step = itmp 263 helium_env(k)%helium%num_steps = helium_env(k)%helium%last_step & 264 - helium_env(k)%helium%first_step 265 ELSE 266 CALL section_vals_val_get(input, "MOTION%PINT%NUM_STEPS", & 267 i_val=itmp) 268 helium_env(k)%helium%num_steps = itmp 269 helium_env(k)%helium%last_step = helium_env(k)%helium%first_step & 270 + helium_env(k)%helium%num_steps 271 END IF 272 273 ! boundary conditions 274 CALL section_vals_val_get(helium_section, "PERIODIC", & 275 l_val=helium_env(k)%helium%periodic) 276 CALL section_vals_val_get(helium_section, "CELL_SHAPE", & 277 i_val=helium_env(k)%helium%cell_shape) 278 279 CALL section_vals_val_get(helium_section, "DROPLET_RADIUS", & 280 r_val=helium_env(k)%helium%droplet_radius) 281 282 ! Set density Rho, number of atoms N and volume V ( Rho = N / V ). 283 ! Allow only 2 out of 3 values to be defined at the same time, calculate 284 ! the third. 285 ! Note, that DENSITY and NATOMS keywords have default values, while 286 ! CELL_SIZE does not. Thus if CELL_SIZE is given explicitly then one and 287 ! only one of the two remaining options must be give explicitly as well. 288 ! If CELL_SIZE is not given explicitly then all four combinations of the 289 ! two other options are valid. 290 CALL section_vals_val_get(helium_section, "DENSITY", & 291 explicit=expl_dens, r_val=helium_env(k)%helium%density) 292 CALL section_vals_val_get(helium_section, "NATOMS", & 293 explicit=expl_nats, i_val=helium_env(k)%helium%atoms) 294 CALL section_vals_val_get(helium_section, "CELL_SIZE", & 295 explicit=expl_cell) 296 cgeof = 1.0_dp 297 IF (helium_env(k)%helium%periodic) THEN 298 IF (helium_env(k)%helium%cell_shape .EQ. helium_cell_shape_octahedron) cgeof = 2.0_dp 299 END IF 300 rtmp = (cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%density)**(1.0_dp/3.0_dp) 301 IF (.NOT. expl_cell) THEN 302 helium_env(k)%helium%cell_size = rtmp 303 ELSE 304 CALL section_vals_val_get(helium_section, "CELL_SIZE", & 305 r_val=helium_env(k)%helium%cell_size) 306 ! only more work if not all three values are consistent: 307 IF (ABS(helium_env(k)%helium%cell_size - rtmp) .GT. 100.0_dp*EPSILON(0.0_dp)* & 308 (ABS(helium_env(k)%helium%cell_size) + rtmp)) THEN 309 IF (expl_dens .AND. expl_nats) THEN 310 msg_str = "DENSITY, NATOMS and CELL_SIZE options "// & 311 "contradict each other" 312 CPWARN(msg_str) 313 END IF 314 !ok we have enough freedom to resolve the conflict: 315 IF (.NOT. expl_dens) THEN 316 helium_env(k)%helium%density = cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%cell_size**3.0_dp 317 IF (.NOT. expl_nats) THEN 318 msg_str = "CELL_SIZE defined but neither "// & 319 "NATOMS nor DENSITY given, using default NATOMS." 320 CPWARN(msg_str) 321 END IF 322 ELSE ! ( expl_dens .AND. .NOT. expl_nats ) 323 ! calculate the nearest number of atoms for given conditions 324 helium_env(k)%helium%atoms = NINT(helium_env(k)%helium%density* & 325 helium_env(k)%helium%cell_size**3.0_dp/cgeof) 326 ! adjust cell size to maintain correct density 327 ! (should be a small correction) 328 rtmp = (cgeof*helium_env(k)%helium%atoms/helium_env(k)%helium%density & 329 )**(1.0_dp/3.0_dp) 330 IF (ABS(helium_env(k)%helium%cell_size - rtmp) .GT. 100.0_dp*EPSILON(0.0_dp) & 331 *(ABS(helium_env(k)%helium%cell_size) + rtmp)) THEN 332 msg_str = "Adjusting actual cell size "// & 333 "to maintain correct density." 334 CPWARN(msg_str) 335 helium_env(k)%helium%cell_size = rtmp 336 END IF 337 END IF 338 END IF 339 END IF 340 helium_env(k)%helium%cell_size_inv = 1.0_dp/helium_env(k)%helium%cell_size 341 ! From now on helium%density, helium%atoms and helium%cell_size are 342 ! correctly defined. 343 344 ! set the M matrix for winding number calculations 345 SELECT CASE (helium_env(k)%helium%cell_shape) 346 347 CASE (helium_cell_shape_octahedron) 348 helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size 349 helium_env(k)%helium%cell_m(2, 1) = 0.0_dp 350 helium_env(k)%helium%cell_m(3, 1) = 0.0_dp 351 helium_env(k)%helium%cell_m(1, 2) = 0.0_dp 352 helium_env(k)%helium%cell_m(2, 2) = helium_env(k)%helium%cell_size 353 helium_env(k)%helium%cell_m(3, 2) = 0.0_dp 354 helium_env(k)%helium%cell_m(1, 3) = helium_env(k)%helium%cell_size/2.0_dp 355 helium_env(k)%helium%cell_m(2, 3) = helium_env(k)%helium%cell_size/2.0_dp 356 helium_env(k)%helium%cell_m(3, 3) = helium_env(k)%helium%cell_size/2.0_dp 357 helium_env(k)%helium%cell_m_inv(1, 1) = 1.0_dp/helium_env(k)%helium%cell_size 358 helium_env(k)%helium%cell_m_inv(2, 1) = 0.0_dp 359 helium_env(k)%helium%cell_m_inv(3, 1) = 0.0_dp 360 helium_env(k)%helium%cell_m_inv(1, 2) = 0.0_dp 361 helium_env(k)%helium%cell_m_inv(2, 2) = 1.0_dp/helium_env(k)%helium%cell_size 362 helium_env(k)%helium%cell_m_inv(3, 2) = 0.0_dp 363 helium_env(k)%helium%cell_m_inv(1, 3) = -1.0_dp/helium_env(k)%helium%cell_size 364 helium_env(k)%helium%cell_m_inv(2, 3) = -1.0_dp/helium_env(k)%helium%cell_size 365 helium_env(k)%helium%cell_m_inv(3, 3) = 2.0_dp/helium_env(k)%helium%cell_size 366 CASE (helium_cell_shape_cube) 367 368 helium_env(k)%helium%cell_m(1, 1) = helium_env(k)%helium%cell_size 369 helium_env(k)%helium%cell_m(2, 1) = 0.0_dp 370 helium_env(k)%helium%cell_m(3, 1) = 0.0_dp 371 helium_env(k)%helium%cell_m(1, 2) = 0.0_dp 372 helium_env(k)%helium%cell_m(2, 2) = helium_env(k)%helium%cell_size 373 helium_env(k)%helium%cell_m(3, 2) = 0.0_dp 374 helium_env(k)%helium%cell_m(1, 3) = 0.0_dp 375 helium_env(k)%helium%cell_m(2, 3) = 0.0_dp 376 helium_env(k)%helium%cell_m(3, 3) = helium_env(k)%helium%cell_size 377 helium_env(k)%helium%cell_m_inv(1, 1) = 1.0_dp/helium_env(k)%helium%cell_size 378 helium_env(k)%helium%cell_m_inv(2, 1) = 0.0_dp 379 helium_env(k)%helium%cell_m_inv(3, 1) = 0.0_dp 380 helium_env(k)%helium%cell_m_inv(1, 2) = 0.0_dp 381 helium_env(k)%helium%cell_m_inv(2, 2) = 1.0_dp/helium_env(k)%helium%cell_size 382 helium_env(k)%helium%cell_m_inv(3, 2) = 0.0_dp 383 helium_env(k)%helium%cell_m_inv(1, 3) = 0.0_dp 384 helium_env(k)%helium%cell_m_inv(2, 3) = 0.0_dp 385 helium_env(k)%helium%cell_m_inv(3, 3) = 1.0_dp/helium_env(k)%helium%cell_size 386 CASE DEFAULT 387 helium_env(k)%helium%cell_m(:, :) = 0.0_dp 388 helium_env(k)%helium%cell_m_inv(:, :) = 0.0_dp 389 390 END SELECT 391 392 END DO ! k 393 394 IF (logger%para_env%ionode) THEN 395 CALL section_vals_val_get(helium_section, "POTENTIAL_FILE_NAME", & 396 c_val=potential_file_name) 397 CALL open_file(file_name=TRIM(potential_file_name), & 398 file_action="READ", file_status="OLD", unit_number=input_unit) 399 READ (input_unit, "(A)") msg_str 400 READ (msg_str, *, IOSTAT=i) nlines, pdx, tau, & 401 x_spline, dx, he_mass 402 IF (i /= 0) THEN 403 ! old style potential file, use default mass and potential 404 he_mass = 4.00263037059764_dp !< 4He mass in [u] 405 expl_pot = .FALSE. 406 READ (msg_str, *, IOSTAT=i) nlines, pdx, tau, & 407 x_spline, dx 408 IF (i /= 0) THEN 409 msg_str = "Format/Read Error from Solvent POTENTIAL_FILE" 410 CPABORT(msg_str) 411 END IF 412 ELSE 413 expl_pot = .TRUE. 414 ! in file really hb2m in kelvin angstrom**2 415 he_mass = angstrom**2*kelvin/massunit/he_mass 416 ! tau might be negative to get older versions of cp2k, 417 ! that cannot handle the new potential file format to 418 ! crash and not run a calculation with wrong mass/potential 419 tau = ABS(tau) 420 END IF 421 tau = kelvin/tau 422 x_spline = x_spline/angstrom 423 dx = dx/angstrom 424 END IF 425 426 CALL mp_bcast(nlines, logger%para_env%source, & 427 helium_env(1)%comm) 428 CALL mp_bcast(pdx, logger%para_env%source, & 429 helium_env(1)%comm) 430 CALL mp_bcast(tau, logger%para_env%source, & 431 helium_env(1)%comm) 432 CALL mp_bcast(x_spline, logger%para_env%source, & 433 helium_env(1)%comm) 434 CALL mp_bcast(dx, logger%para_env%source, & 435 helium_env(1)%comm) 436 CALL mp_bcast(he_mass, logger%para_env%source, & 437 helium_env(1)%comm) 438 isize = (pdx + 1)*(pdx + 2) + 1 439 ALLOCATE (pot_transfer(nlines, isize)) 440 IF (logger%para_env%ionode) THEN 441 IF (expl_pot) THEN 442 DO i = 1, nlines 443 READ (input_unit, *) pot_transfer(i, :) 444 END DO 445 ELSE 446 DO i = 1, nlines 447 READ (input_unit, *) pot_transfer(i, 1:isize - 1) 448 ! potential implicit, calculate it here now: 449 pot_transfer(i, isize) = helium_vij(x_spline + (i - 1)*dx)*kelvin 450 END DO 451 END IF 452 CALL close_file(unit_number=input_unit) 453 END IF 454 CALL mp_bcast(pot_transfer, logger%para_env%source, & 455 helium_env(1)%comm) 456 457 CALL spline_data_create(helium_env(1)%helium%vij) 458 CALL init_splinexy(helium_env(1)%helium%vij, nlines) 459 helium_env(1)%helium%vij%x1 = x_spline 460 461 CALL spline_data_create(helium_env(1)%helium%u0) 462 CALL init_splinexy(helium_env(1)%helium%u0, nlines) 463 helium_env(1)%helium%u0%x1 = x_spline 464 465 CALL spline_data_create(helium_env(1)%helium%e0) 466 CALL init_splinexy(helium_env(1)%helium%e0, nlines) 467 helium_env(1)%helium%e0%x1 = x_spline 468 469 isize = pdx + 1 470 ntab = ((isize + 1)*isize)/2 - 1 ! -1 because endpoint approx treated separately 471 ALLOCATE (helium_env(1)%helium%uoffdiag(ntab, 2, nlines)) 472 ALLOCATE (helium_env(1)%helium%eoffdiag(ntab, 2, nlines)) 473 DO j = 1, isize 474 DO i = j, isize 475 IF (i + j == 2) CYCLE ! endpoint approx later separately 476 k = ((i - 1)*i)/2 + j 477 helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2) 478 CALL init_spline(helium_env(1)%helium%vij, dx=dx) 479 helium_env(1)%helium%uoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:) 480 helium_env(1)%helium%uoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:) 481 k = k + ((isize + 1)*isize)/2 482 helium_env(1)%helium%vij%y(:) = pot_transfer(:, k)*angstrom**(2*i - 2)/kelvin 483 CALL init_spline(helium_env(1)%helium%vij, dx=dx) 484 helium_env(1)%helium%eoffdiag(ntab, 1, :) = helium_env(1)%helium%vij%y(:) 485 helium_env(1)%helium%eoffdiag(ntab, 2, :) = helium_env(1)%helium%vij%y2(:) 486 ntab = ntab - 1 487 END DO 488 END DO 489 490 ntab = SIZE(pot_transfer, 2) 491 helium_env(1)%helium%vij%y(:) = pot_transfer(:, ntab)/kelvin 492 CALL init_spline(helium_env(1)%helium%vij, dx=dx) 493 494 helium_env(1)%helium%u0%y(:) = pot_transfer(:, 1) 495 CALL init_spline(helium_env(1)%helium%u0, dx=dx) 496 k = ((isize + 1)*isize)/2 + 1 497 helium_env(1)%helium%e0%y(:) = pot_transfer(:, k)/kelvin 498 CALL init_spline(helium_env(1)%helium%e0, dx=dx) 499 500 DO k = 2, mepos 501 helium_env(k)%helium%vij => helium_env(1)%helium%vij 502 helium_env(k)%helium%u0 => helium_env(1)%helium%u0 503 helium_env(k)%helium%e0 => helium_env(1)%helium%e0 504 helium_env(k)%helium%uoffdiag => helium_env(1)%helium%uoffdiag 505 helium_env(k)%helium%eoffdiag => helium_env(1)%helium%eoffdiag 506 END DO 507 508 DO k = 1, mepos 509 510 helium_env(k)%helium%pdx = pdx 511 helium_env(k)%helium%tau = tau 512 513 ! boltzmann : Boltzmann constant [J/K] 514 ! h_bar : Planck constant [J*s] 515 ! J = kg*m^2/s^2 516 ! 4He mass in [kg] 517 mHe = he_mass*a_mass 518 ! physical temperature [K] 519 T = kelvin/helium_env(k)%helium%tau/helium_env(k)%helium%beads 520 ! prefactors for calculating superfluid fractions [Angstrom^-2] 521 helium_env(k)%helium%wpref = (((1e-20/h_bar)*mHe)/h_bar)*boltzmann*T 522 helium_env(k)%helium%apref = (((4e-20/h_bar)*mHe)/h_bar)*boltzmann*T 523 524 helium_env(k)%helium%he_mass_au = he_mass*massunit 525 helium_env(k)%helium%hb2m = 1.0_dp/helium_env(k)%helium%he_mass_au 526 helium_env(k)%helium%pweight = 0.0_dp 527 528 ! Choose sampling method: 529 CALL section_vals_val_get(helium_section, "SAMPLING_METHOD", & 530 i_val=helium_env(k)%helium%sampling_method) 531 532 SELECT CASE (helium_env(k)%helium%sampling_method) 533 CASE (helium_sampling_ceperley) 534 ! check value of maxcycle 535 CALL section_vals_val_get(helium_section, "CEPERLEY%MAX_PERM_CYCLE", & 536 i_val=helium_env(k)%helium%maxcycle) 537 i = helium_env(k)%helium%maxcycle 538 CPASSERT(i >= 0) 539 i = helium_env(k)%helium%atoms - helium_env(k)%helium%maxcycle 540 CPASSERT(i >= 0) 541 542 ! set m-distribution parameters 543 CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%DISTRIBUTION-TYPE", & 544 i_val=i) 545 CPASSERT(i >= 1) 546 CPASSERT(i <= 6) 547 helium_env(k)%helium%m_dist_type = i 548 CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-VALUE", & 549 i_val=i) 550 CPASSERT(i >= 1) 551 CPASSERT(i <= helium_env(k)%helium%maxcycle) 552 helium_env(k)%helium%m_value = i 553 CALL section_vals_val_get(helium_section, "CEPERLEY%M-SAMPLING%M-RATIO", & 554 r_val=rtmp) 555 CPASSERT(rtmp > 0.0_dp) 556 CPASSERT(rtmp <= 1.0_dp) 557 helium_env(k)%helium%m_ratio = rtmp 558 559 CALL section_vals_val_get(helium_section, "CEPERLEY%BISECTION", & 560 i_val=helium_env(k)%helium%bisection) 561 ! precheck bisection value (not all invalids are filtered out here yet) 562 i = helium_env(k)%helium%bisection 563 CPASSERT(i > 1) 564 i = helium_env(k)%helium%beads - helium_env(k)%helium%bisection 565 CPASSERT(i > 0) 566 ! 567 itmp = helium_env(k)%helium%bisection 568 rtmp = 2.0_dp**(ANINT(LOG(REAL(itmp, dp))/LOG(2.0_dp))) 569 tcheck = ABS(REAL(itmp, KIND=dp) - rtmp) 570 IF (tcheck > 100.0_dp*EPSILON(0.0_dp)) THEN 571 msg_str = "BISECTION should be integer power of 2." 572 CPABORT(msg_str) 573 END IF 574 helium_env(k)%helium%bisctlog2 = NINT(LOG(REAL(itmp, dp))/LOG(2.0_dp)) 575 576 CASE (helium_sampling_worm) 577 NULLIFY (input_worm) 578 input_worm => section_vals_get_subs_vals(helium_env(k)%helium%input, & 579 "MOTION%PINT%HELIUM%WORM") 580 CALL section_vals_val_get(helium_section, "WORM%CENTROID_DRMAX", & 581 r_val=helium_env(k)%helium%worm_centroid_drmax) 582 583 CALL section_vals_val_get(helium_section, "WORM%STAGING_L", & 584 i_val=helium_env(k)%helium%worm_staging_l) 585 586 CALL section_vals_val_get(helium_section, "WORM%OPEN_CLOSE_SCALE", & 587 r_val=helium_env(k)%helium%worm_open_close_scale) 588 589 CALL section_vals_val_get(helium_section, "WORM%ALLOW_OPEN", & 590 l_val=helium_env(k)%helium%worm_allow_open) 591 592 IF (helium_env(k)%helium%worm_staging_l + 1 >= helium_env(k)%helium%beads) THEN 593 msg_str = "STAGING_L for worm sampling is to large" 594 CPABORT(msg_str) 595 ELSE IF (helium_env(k)%helium%worm_staging_l < 1) THEN 596 msg_str = "STAGING_L must be positive integer" 597 CPABORT(msg_str) 598 END IF 599 600 CALL section_vals_val_get(helium_section, "WORM%SHOW_STATISTICS", & 601 l_val=helium_env(k)%helium%worm_show_statistics) 602 603 ! precompute an expensive scaling for the open and close acceptance probability 604 ! tau is not included here, as It will be first defined in a few lines 605 rtmp = 2.0_dp*pi*helium_env(k)%helium%hb2m 606 rtmp = SQRT(rtmp) 607 rtmp = rtmp**3 608 rtmp = rtmp*helium_env(k)%helium%worm_open_close_scale 609 IF (helium_env(k)%helium%periodic) THEN 610 rtmp = rtmp*helium_env(k)%helium%density 611 ELSE 612 rtmp = rtmp*helium_env(k)%helium%atoms/ & 613 (4.0_dp/3.0_dp*pi*helium_env(k)%helium%droplet_radius**3) 614 END IF 615 helium_env(k)%helium%worm_ln_openclose_scale = LOG(rtmp) 616 617 ! deal with acceptance statistics without changing the ceperley stuff 618 helium_env(k)%helium%maxcycle = 1 619 helium_env(k)%helium%bisctlog2 = 0 620 621 ! get the absolute weights of the individual moves 622 helium_env(k)%helium%worm_all_limit = 0 623 CALL section_vals_val_get(helium_section, "WORM%CENTROID_WEIGHT", & 624 i_val=itmp) 625 helium_env(k)%helium%worm_centroid_min = 1 626 helium_env(k)%helium%worm_centroid_max = itmp 627 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 628 CALL section_vals_val_get(helium_section, "WORM%STAGING_WEIGHT", & 629 i_val=itmp) 630 helium_env(k)%helium%worm_staging_min = helium_env(k)%helium%worm_centroid_max + 1 631 helium_env(k)%helium%worm_staging_max = helium_env(k)%helium%worm_centroid_max + itmp 632 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 633 IF (helium_env(k)%helium%worm_allow_open) THEN 634 CALL section_vals_val_get(helium_section, "WORM%CRAWL_WEIGHT", & 635 i_val=itmp) 636 helium_env(k)%helium%worm_fcrawl_min = helium_env(k)%helium%worm_staging_max + 1 637 helium_env(k)%helium%worm_fcrawl_max = helium_env(k)%helium%worm_staging_max + itmp 638 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 639 helium_env(k)%helium%worm_bcrawl_min = helium_env(k)%helium%worm_fcrawl_max + 1 640 helium_env(k)%helium%worm_bcrawl_max = helium_env(k)%helium%worm_fcrawl_max + itmp 641 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 642 CALL section_vals_val_get(helium_section, "WORM%HEAD_TAIL_WEIGHT", & 643 i_val=itmp) 644 helium_env(k)%helium%worm_head_min = helium_env(k)%helium%worm_bcrawl_max + 1 645 helium_env(k)%helium%worm_head_max = helium_env(k)%helium%worm_bcrawl_max + itmp 646 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 647 helium_env(k)%helium%worm_tail_min = helium_env(k)%helium%worm_head_max + 1 648 helium_env(k)%helium%worm_tail_max = helium_env(k)%helium%worm_head_max + itmp 649 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 650 CALL section_vals_val_get(helium_section, "WORM%SWAP_WEIGHT", & 651 i_val=itmp) 652 helium_env(k)%helium%worm_swap_min = helium_env(k)%helium%worm_tail_max + 1 653 helium_env(k)%helium%worm_swap_max = helium_env(k)%helium%worm_tail_max + itmp 654 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 655 CALL section_vals_val_get(helium_section, "WORM%OPEN_CLOSE_WEIGHT", & 656 i_val=itmp) 657 helium_env(k)%helium%worm_open_close_min = helium_env(k)%helium%worm_swap_max + 1 658 helium_env(k)%helium%worm_open_close_max = helium_env(k)%helium%worm_swap_max + itmp 659 helium_env(k)%helium%worm_all_limit = helium_env(k)%helium%worm_all_limit + itmp 660 CALL section_vals_val_get(helium_section, "WORM%CRAWL_REPETITION", & 661 i_val=helium_env(k)%helium%worm_repeat_crawl) 662 END IF 663 664 !CPPostcondition(i<helium_env(k)%helium%beads,cp_failure_level,routineP,failure) 665 ! end of worm 666 CASE DEFAULT 667 msg_str = "Unknown helium sampling method!" 668 CPABORT(msg_str) 669 END SELECT 670 671 ! ALLOCATE helium-related arrays 672 i = helium_env(k)%helium%atoms 673 j = helium_env(k)%helium%beads 674 ALLOCATE (helium_env(k)%helium%pos(3, i, j)) 675 helium_env(k)%helium%pos = 0.0_dp 676 ALLOCATE (helium_env(k)%helium%work(3, i, j)) 677 ALLOCATE (helium_env(k)%helium%ptable(helium_env(k)%helium%maxcycle + 1)) 678 ALLOCATE (helium_env(k)%helium%permutation(i)) 679 ALLOCATE (helium_env(k)%helium%iperm(i)) 680 ALLOCATE (helium_env(k)%helium%tmatrix(i, i)) 681 ALLOCATE (helium_env(k)%helium%nmatrix(i, 2*i)) 682 ALLOCATE (helium_env(k)%helium%pmatrix(i, i)) 683 ALLOCATE (helium_env(k)%helium%ipmatrix(i, i)) 684 itmp = helium_env(k)%helium%bisctlog2 + 2 685 ALLOCATE (helium_env(k)%helium%num_accepted(itmp, helium_env(k)%helium%maxcycle)) 686 ALLOCATE (helium_env(k)%helium%plength_avrg(helium_env(k)%helium%atoms)) 687 ALLOCATE (helium_env(k)%helium%plength_inst(helium_env(k)%helium%atoms)) 688 ALLOCATE (helium_env(k)%helium%atom_plength(helium_env(k)%helium%atoms)) 689 690 ! check whether rdfs should be calculated and printed 691 helium_env(k)%helium%rdf_present = helium_property_active(helium_env(k)%helium, "RDF") 692 IF (helium_env(k)%helium%rdf_present) THEN 693 ! allocate & initialize rdf related data structures 694 CALL helium_rdf_init(helium_env(k)%helium) 695 END IF 696 697 ! check whether densities should be calculated and printed 698 helium_env(k)%helium%rho_present = helium_property_active(helium_env(k)%helium, "RHO") 699 IF (helium_env(k)%helium%rho_present) THEN 700 ! allocate & initialize density related data structures 701 NULLIFY (helium_env(k)%helium%rho_property) 702 CALL helium_rho_init(helium_env(k)%helium) 703 END IF 704 705 END DO 706 707 ! restore averages calculated in previous runs 708 CALL helium_averages_restore(helium_env) 709 710 DO k = 1, mepos 711 ! fill in the solute-related data structures 712 helium_env(k)%helium%e_corr = 0.0_dp 713 IF (helium_env(k)%helium%solute_present) THEN 714 IF (helium_env(k)%helium%solute_beads > helium_env(k)%helium%beads) THEN 715 ! Imaginary time striding for solute: 716 helium_env(k)%helium%bead_ratio = helium_env(k)%helium%solute_beads/ & 717 helium_env(k)%helium%beads 718 ! check if bead numbers are commensurate: 719 i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%beads - helium_env(k)%helium%solute_beads 720 IF (i /= 0) THEN 721 msg_str = "Adjust number of solute beads to multiple of solvent beads." 722 CPABORT(msg_str) 723 END IF 724 msg_str = "Using multiple-time stepping in imaginary time for solute to couple "// & 725 "to fewer solvent beads, e.g. for factor 3: "// & 726 "he_1 - 3*sol_1; he_2 - 3*sol_4... "// & 727 "Avoid too large coupling factors." 728 CPWARN(msg_str) 729 ELSE IF (helium_env(k)%helium%solute_beads < helium_env(k)%helium%beads) THEN 730 ! Imaginary time striding for solvent: 731 helium_env(k)%helium%bead_ratio = helium_env(k)%helium%beads/ & 732 helium_env(k)%helium%solute_beads 733 ! check if bead numbers are commensurate: 734 i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%solute_beads - helium_env(k)%helium%beads 735 IF (i /= 0) THEN 736 msg_str = "Adjust number of solvent beads to multiple of solute beads." 737 CPABORT(msg_str) 738 END IF 739 msg_str = "Coupling solvent beads to fewer solute beads via "// & 740 "direct coupling, e.g. for factor 3: "// & 741 "sol_1 - he_1,2,3; sol_2 - he_4,5,6..." 742 CPWARN(msg_str) 743 END IF 744!TODO Adjust helium bead number if not comm. and if coords not given expl. 745 746 ! check if tau, temperature and bead number are consistent: 747 tcheck = ABS((helium_env(k)%helium%tau*helium_env(k)%helium%beads - solute%beta)/solute%beta) 748 IF (tcheck > 1.0e-14_dp) THEN 749 msg_str = "Tau, temperature and bead number are inconsistent." 750 CPABORT(msg_str) 751 END IF 752 753 CALL helium_set_solute_indices(helium_env(k)%helium, solute) 754 CALL helium_set_solute_cell(helium_env(k)%helium, solute) 755 756 ! set the interaction potential type 757 CALL section_vals_val_get(helium_section, "SOLUTE_INTERACTION", & 758 i_val=helium_env(k)%helium%solute_interaction) 759 IF (helium_env(k)%helium%solute_interaction .EQ. helium_solute_intpot_none) THEN 760 WRITE (msg_str, '(A,I0,A)') & 761 "Solute found but no helium-solute interaction selected "// & 762 "(see SOLUTE_INTERACTION keyword)" 763 CPABORT(msg_str) 764 END IF 765 766 ! ALLOCATE solute-related arrays 767 ALLOCATE (helium_env(k)%helium%force_avrg(helium_env(k)%helium%solute_beads, & 768 helium_env(k)%helium%solute_atoms*3)) 769 ALLOCATE (helium_env(k)%helium%force_inst(helium_env(k)%helium%solute_beads, & 770 helium_env(k)%helium%solute_atoms*3)) 771 772 ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_1d(solute%p*solute%ndim)) 773 ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_np_1d(solute%p*solute%ndim*helium_env(k)%helium%num_env)) 774 ALLOCATE (helium_env(k)%helium%rtmp_p_ndim_2d(solute%p, solute%ndim)) 775 776 ELSE 777 helium_env(k)%helium%bead_ratio = 0 778 IF (helium_env(k)%helium%periodic) THEN 779 ! this assumes a specific potential (and its ugly): 780 x1 = angstrom*0.5_dp*helium_env(k)%helium%cell_size 781 ! 10.8 is in Kelvin, x1 needs to be in Angstrom, 782 ! since 2.9673 is in Angstrom 783 helium_env(k)%helium%e_corr = (twopi*helium_env(k)%helium%density/angstrom**3*10.8_dp* & 784 (544850.4_dp*EXP(-13.353384_dp*x1/2.9673_dp)*(2.9673_dp/13.353384_dp)**3* & 785 (2.0_dp + 2.0_dp*13.353384_dp*x1/2.9673_dp + (13.353384_dp*x1/2.9673_dp)**2) - & 786 (((0.1781_dp/7.0_dp*(2.9673_dp/x1)**2 + 0.4253785_dp/5.0_dp)*(2.9673_dp/x1)**2 + & 787 1.3732412_dp/3.0_dp)*(2.9673_dp/x1)**3)*2.9673_dp**3))/kelvin 788 END IF 789 END IF 790 791 ! ALLOCATE temporary arrays 792 ALLOCATE (helium_env(k)%helium%itmp_atoms_1d(helium_env(k)%helium%atoms)) 793 ALLOCATE (helium_env(k)%helium%ltmp_atoms_1d(helium_env(k)%helium%atoms)) 794 ALLOCATE (helium_env(k)%helium%itmp_atoms_np_1d(helium_env(k)%helium%atoms*helium_env(k)%helium%num_env)) 795 ALLOCATE (helium_env(k)%helium%rtmp_3_np_1d(3*helium_env(k)%helium%num_env)) 796 ALLOCATE (helium_env(k)%helium%rtmp_3_atoms_beads_1d(3*helium_env(k)%helium%atoms* & 797 helium_env(k)%helium%beads)) 798 ALLOCATE (helium_env(k)%helium%rtmp_3_atoms_beads_np_1d(3*helium_env(k)%helium%atoms* & 799 helium_env(k)%helium%beads*helium_env(k)%helium%num_env)) 800 ALLOCATE (helium_env(k)%helium%ltmp_3_atoms_beads_3d(3, helium_env(k)%helium%atoms, & 801 helium_env(k)%helium%beads)) 802 IF (k .EQ. 1) THEN 803 CALL helium_write_setup(helium_env(k)%helium) 804 END IF 805 END DO 806 DEALLOCATE (pot_transfer) 807 ELSE 808 ! Deallocate env_all on processors without helium_env 809 DEALLOCATE (env_all) 810 END IF ! mepos .GT. 0 811 812 NULLIFY (env_all) 813 CALL timestop(handle) 814 815 RETURN 816 END SUBROUTINE helium_create 817 818! *************************************************************************** 819!> \brief Releases helium_solvent_type 820!> \param helium_env ... 821!> \par History 822!> 2016-07-14 Modified to work with independent helium_env [cschran] 823!> \author hforbert 824! ************************************************************************************************** 825 SUBROUTINE helium_release(helium_env) 826 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 827 828 CHARACTER(len=*), PARAMETER :: routineN = 'helium_release', routineP = moduleN//':'//routineN 829 830 INTEGER :: k 831 832 IF (ASSOCIATED(helium_env)) THEN 833 DO k = 1, SIZE(helium_env) 834 IF (k .EQ. 1) THEN 835 CALL mp_comm_free(helium_env(k)%comm) 836 DEALLOCATE (helium_env(k)%env_all) 837 END IF 838 NULLIFY (helium_env(k)%env_all) 839 helium_env(k)%helium%ref_count = helium_env(k)%helium%ref_count - 1 840 IF (helium_env(k)%helium%ref_count < 1) THEN 841 842 ! DEALLOCATE temporary arrays 843 DEALLOCATE ( & 844 helium_env(k)%helium%ltmp_3_atoms_beads_3d, & 845 helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, & 846 helium_env(k)%helium%rtmp_3_atoms_beads_1d, & 847 helium_env(k)%helium%rtmp_3_np_1d, & 848 helium_env(k)%helium%itmp_atoms_np_1d, & 849 helium_env(k)%helium%ltmp_atoms_1d, & 850 helium_env(k)%helium%itmp_atoms_1d) 851 852 NULLIFY ( & 853 helium_env(k)%helium%ltmp_3_atoms_beads_3d, & 854 helium_env(k)%helium%rtmp_3_atoms_beads_np_1d, & 855 helium_env(k)%helium%rtmp_3_atoms_beads_1d, & 856 helium_env(k)%helium%rtmp_3_np_1d, & 857 helium_env(k)%helium%itmp_atoms_np_1d, & 858 helium_env(k)%helium%ltmp_atoms_1d, & 859 helium_env(k)%helium%itmp_atoms_1d & 860 ) 861 862 IF (helium_env(k)%helium%solute_present) THEN 863 ! DEALLOCATE solute-related arrays 864 DEALLOCATE ( & 865 helium_env(k)%helium%rtmp_p_ndim_2d, & 866 helium_env(k)%helium%rtmp_p_ndim_np_1d, & 867 helium_env(k)%helium%rtmp_p_ndim_1d, & 868 helium_env(k)%helium%force_inst, & 869 helium_env(k)%helium%force_avrg) 870 NULLIFY ( & 871 helium_env(k)%helium%rtmp_p_ndim_2d, & 872 helium_env(k)%helium%rtmp_p_ndim_np_1d, & 873 helium_env(k)%helium%rtmp_p_ndim_1d, & 874 helium_env(k)%helium%force_inst, & 875 helium_env(k)%helium%force_avrg) 876 END IF 877 878 IF (helium_env(k)%helium%rho_present) THEN 879 DEALLOCATE ( & 880 helium_env(k)%helium%rho_rstr, & 881 helium_env(k)%helium%rho_accu, & 882 helium_env(k)%helium%rho_inst, & 883 helium_env(k)%helium%rho_incr) 884 NULLIFY ( & 885 helium_env(k)%helium%rho_rstr, & 886 helium_env(k)%helium%rho_accu, & 887 helium_env(k)%helium%rho_inst, & 888 helium_env(k)%helium%rho_incr) 889 ! DEALLOCATE everything 890 DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix) 891 DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_name) 892 DEALLOCATE (helium_env(k)%helium%rho_property(rho_atom_number)%component_index) 893 NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%filename_suffix) 894 NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_name) 895 NULLIFY (helium_env(k)%helium%rho_property(rho_atom_number)%component_index) 896 DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix) 897 DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_name) 898 DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_number)%component_index) 899 NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%filename_suffix) 900 NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_name) 901 NULLIFY (helium_env(k)%helium%rho_property(rho_winding_number)%component_index) 902 DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix) 903 DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name) 904 DEALLOCATE (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index) 905 NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%filename_suffix) 906 NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_name) 907 NULLIFY (helium_env(k)%helium%rho_property(rho_winding_cycle)%component_index) 908 DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix) 909 DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_name) 910 DEALLOCATE (helium_env(k)%helium%rho_property(rho_projected_area)%component_index) 911 NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%filename_suffix) 912 NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_name) 913 NULLIFY (helium_env(k)%helium%rho_property(rho_projected_area)%component_index) 914 DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix) 915 DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name) 916 DEALLOCATE (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index) 917 NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%filename_suffix) 918 NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_name) 919 NULLIFY (helium_env(k)%helium%rho_property(rho_moment_of_inertia)%component_index) 920 DEALLOCATE (helium_env(k)%helium%rho_property) 921 NULLIFY (helium_env(k)%helium%rho_property) 922 END IF 923 924 CALL helium_rdf_release(helium_env(k)%helium) 925 926 ! DEALLOCATE helium-related arrays 927 DEALLOCATE ( & 928 helium_env(k)%helium%atom_plength, & 929 helium_env(k)%helium%plength_inst, & 930 helium_env(k)%helium%plength_avrg, & 931 helium_env(k)%helium%num_accepted, & 932 helium_env(k)%helium%ipmatrix, & 933 helium_env(k)%helium%pmatrix, & 934 helium_env(k)%helium%nmatrix, & 935 helium_env(k)%helium%tmatrix, & 936 helium_env(k)%helium%iperm, & 937 helium_env(k)%helium%permutation, & 938 helium_env(k)%helium%ptable, & 939 helium_env(k)%helium%work, & 940 helium_env(k)%helium%pos) 941 NULLIFY ( & 942 helium_env(k)%helium%atom_plength, & 943 helium_env(k)%helium%plength_inst, & 944 helium_env(k)%helium%plength_avrg, & 945 helium_env(k)%helium%num_accepted, & 946 helium_env(k)%helium%ipmatrix, & 947 helium_env(k)%helium%pmatrix, & 948 helium_env(k)%helium%nmatrix, & 949 helium_env(k)%helium%tmatrix, & 950 helium_env(k)%helium%iperm, & 951 helium_env(k)%helium%permutation, & 952 helium_env(k)%helium%ptable, & 953 helium_env(k)%helium%work, & 954 helium_env(k)%helium%pos & 955 ) 956 957 IF (k == 1) THEN 958 CALL spline_data_release(helium_env(k)%helium%vij) 959 CALL spline_data_release(helium_env(k)%helium%u0) 960 CALL spline_data_release(helium_env(k)%helium%e0) 961 DEALLOCATE (helium_env(k)%helium%uoffdiag) 962 DEALLOCATE (helium_env(k)%helium%eoffdiag) 963 END IF 964 NULLIFY (helium_env(k)%helium%uoffdiag, & 965 helium_env(k)%helium%eoffdiag, & 966 helium_env(k)%helium%vij, & 967 helium_env(k)%helium%u0, & 968 helium_env(k)%helium%e0) 969 970 CALL delete_rng_stream(helium_env(k)%helium%rng_stream_uniform) 971 CALL delete_rng_stream(helium_env(k)%helium%rng_stream_gaussian) 972 973 ! deallocate solute-related arrays 974 IF (helium_env(k)%helium%solute_present) THEN 975 DEALLOCATE (helium_env(k)%helium%solute_element) 976 NULLIFY (helium_env(k)%helium%solute_element) 977 END IF 978 979 ! Deallocate everything from the helium_set_solute_indices 980 IF (ASSOCIATED(helium_env(k)%helium%ename)) THEN 981 DEALLOCATE (helium_env(k)%helium%ename) 982 NULLIFY (helium_env(k)%helium%ename) 983 END IF 984 985 DEALLOCATE (helium_env(k)%helium) 986 987 END IF 988 END DO 989 990 DEALLOCATE (helium_env) 991 NULLIFY (helium_env) 992 END IF 993 RETURN 994 END SUBROUTINE helium_release 995 996! *************************************************************************** 997!> \brief Initialize helium data structures. 998!> \param helium_env ... 999!> \param pint_env ... 1000!> \par History 1001!> removed refereces to pint_env_type data structure [lwalewski] 1002!> 2009-11-10 init/restore coords, perm, RNG and forces [lwalewski] 1003!> 2016-07-14 Modified to work with independent helium_env [cschran] 1004!> \author hforbert 1005!> \note Initializes helium coordinates either as random positions or from 1006!> HELIUM%COORD section if it's present in the input file. 1007!> Initializes helium permutation state as identity permutation or 1008!> from HELIUM%PERM section if it's present in the input file. 1009! ************************************************************************************************** 1010 SUBROUTINE helium_init(helium_env, pint_env) 1011 1012 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1013 TYPE(pint_env_type), POINTER :: pint_env 1014 1015 CHARACTER(len=*), PARAMETER :: routineN = 'helium_init', routineP = moduleN//':'//routineN 1016 1017 INTEGER :: handle, k 1018 LOGICAL :: coords_presampled, explicit, presample 1019 REAL(KIND=dp) :: initkT, solute_radius 1020 TYPE(cp_logger_type), POINTER :: logger 1021 TYPE(section_vals_type), POINTER :: helium_section, sec 1022 1023 CALL timeset(routineN, handle) 1024 1025 NULLIFY (logger) 1026 logger => cp_get_default_logger() 1027 1028 IF (ASSOCIATED(helium_env)) THEN 1029 1030 NULLIFY (helium_section) 1031 helium_section => section_vals_get_subs_vals(helium_env(1)%helium%input, & 1032 "MOTION%PINT%HELIUM") 1033 1034 ! restore RNG state 1035 NULLIFY (sec) 1036 sec => section_vals_get_subs_vals(helium_section, "RNG_STATE") 1037 CALL section_vals_get(sec, explicit=explicit) 1038 IF (explicit) THEN 1039 CALL helium_rng_restore(helium_env) 1040 ELSE 1041 CALL helium_write_line("RNG state initialized as new.") 1042 END IF 1043 1044 ! init/restore permutation state 1045 NULLIFY (sec) 1046 sec => section_vals_get_subs_vals(helium_section, "PERM") 1047 CALL section_vals_get(sec, explicit=explicit) 1048 IF (explicit) THEN 1049 CALL helium_perm_restore(helium_env) 1050 ELSE 1051 CALL helium_perm_init(helium_env) 1052 CALL helium_write_line("Permutation state initialized as identity.") 1053 END IF 1054 1055 ! Specify if forces should be obtained as AVG or LAST 1056 DO k = 1, SIZE(helium_env) 1057 CALL section_vals_val_get(helium_section, "GET_FORCES", & 1058 i_val=helium_env(k)%helium%get_helium_forces) 1059 END DO 1060 1061 DO k = 1, SIZE(helium_env) 1062 ! init center of mass 1063 IF (helium_env(k)%helium%solute_present) THEN 1064 helium_env(k)%helium%center(:) = pint_com_pos(pint_env) 1065 ELSE 1066 helium_env(k)%helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/) 1067 END IF 1068 END DO 1069 1070 ! init/restore coordinates 1071 NULLIFY (sec) 1072 sec => section_vals_get_subs_vals(helium_section, "COORD") 1073 CALL section_vals_get(sec, explicit=explicit) 1074 IF (explicit) THEN 1075 CALL helium_coord_restore(helium_env) 1076 CALL helium_write_line("Coordinates restarted.") 1077 ELSE 1078 CALL section_vals_val_get(helium_section, "COORD_INIT_TEMP", r_val=initkT) 1079 CALL section_vals_val_get(helium_section, "SOLUTE_RADIUS", r_val=solute_radius) 1080 CALL helium_coord_init(helium_env, initkT, solute_radius) 1081 IF (initkT > 0.0_dp) THEN 1082 CALL helium_write_line("Coordinates initialized with thermal gaussian.") 1083 ELSE 1084 CALL helium_write_line("Coordinates initialized as point particles.") 1085 END IF 1086 END IF 1087 1088 DO k = 1, SIZE(helium_env) 1089 1090 helium_env(k)%helium%worm_is_closed = .TRUE. 1091 helium_env(k)%helium%worm_atom_idx = 0 1092 helium_env(k)%helium%worm_bead_idx = 0 1093 1094 helium_env(k)%helium%work(:, :, :) = helium_env(k)%helium%pos(:, :, :) 1095 1096 ! init center of mass 1097 IF (helium_env(k)%helium%solute_present) THEN 1098 helium_env(k)%helium%center(:) = pint_com_pos(pint_env) 1099 ELSE 1100 IF (helium_env(k)%helium%periodic) THEN 1101 helium_env(k)%helium%center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/) 1102 ELSE 1103 helium_env(k)%helium%center(:) = helium_com(helium_env(k)%helium) 1104 END IF 1105 END IF 1106 END DO 1107 1108 ! Optional helium coordinate presampling: 1109 ! Assume IONODE to have always at least one helium_env 1110 CALL section_vals_val_get(helium_section, "PRESAMPLE", & 1111 l_val=presample) 1112 coords_presampled = .FALSE. 1113 IF (presample) THEN 1114 DO k = 1, SIZE(helium_env) 1115 helium_env(k)%helium%current_step = 0 1116 END DO 1117 CALL helium_sample(helium_env, pint_env) 1118 DO k = 1, SIZE(helium_env) 1119 IF (helium_env(k)%helium%solute_present) helium_env(k)%helium%force_avrg(:, :) = 0.0_dp 1120 helium_env(k)%helium%energy_avrg(:) = 0.0_dp 1121 helium_env(k)%helium%plength_avrg(:) = 0.0_dp 1122 helium_env(k)%helium%num_accepted(:, :) = 0.0_dp 1123 ! Reset properties accumulated over presample: 1124 helium_env(k)%helium%proarea%accu(:) = 0.0_dp 1125 helium_env(k)%helium%prarea2%accu(:) = 0.0_dp 1126 helium_env(k)%helium%wnmber2%accu(:) = 0.0_dp 1127 helium_env(k)%helium%mominer%accu(:) = 0.0_dp 1128 IF (helium_env(k)%helium%rho_present) THEN 1129 helium_env(k)%helium%rho_accu(:, :, :, :) = 0.0_dp 1130 END IF 1131 IF (helium_env(k)%helium%rdf_present) THEN 1132 helium_env(k)%helium%rdf_accu(:, :) = 0.0_dp 1133 END IF 1134 END DO 1135 coords_presampled = .TRUE. 1136 CALL helium_write_line("Bead coordinates pre-sampled.") 1137 END IF 1138 1139 IF (helium_env(1)%helium%solute_present) THEN 1140 ! restore helium forces 1141 NULLIFY (sec) 1142 sec => section_vals_get_subs_vals(helium_section, "FORCE") 1143 CALL section_vals_get(sec, explicit=explicit) 1144 IF (explicit) THEN 1145 IF (.NOT. coords_presampled) THEN 1146 CALL helium_force_restore(helium_env) 1147 END IF 1148 ELSE 1149 IF (.NOT. coords_presampled) THEN 1150 CALL helium_force_init(helium_env) 1151 CALL helium_write_line("Forces on the solute initialized as zero.") 1152 END IF 1153 END IF 1154 !! Update pint_env force, assume always one helium_env at IONODE 1155 !IF (pint_env%logger%para_env%ionode) THEN 1156 ! pint_env%f(:, :) = pint_env%f(:, :) + helium_env(1)%helium%force_avrg(:, :) 1157 !END IF 1158 !CALL mp_bcast(pint_env%f,& 1159 ! pint_env%logger%para_env%source,& 1160 ! pint_env%logger%para_env%group) 1161 1162 END IF 1163 END IF 1164 1165 CALL timestop(handle) 1166 1167 RETURN 1168 END SUBROUTINE helium_init 1169 1170! *************************************************************************** 1171! Data transfer functions. 1172! 1173! These functions manipulate and transfer data between the runtime 1174! environment and the input structure. 1175! *************************************************************************** 1176 1177! *************************************************************************** 1178!> \brief Initialize helium coordinates with random positions. 1179!> \param helium_env ... 1180!> \param initkT ... 1181!> \param solute_radius ... 1182!> \date 2009-11-09 1183!> \par History 1184!> 2016-07-14 Modified to work with independent helium_env [cschran] 1185!> 2018-04-30 Usefull initialization for droplets [fuhl] 1186!> \author Lukasz Walewski 1187! ************************************************************************************************** 1188 SUBROUTINE helium_coord_init(helium_env, initkT, solute_radius) 1189 TYPE(helium_solvent_p_type), DIMENSION(:), & 1190 INTENT(INOUT), POINTER :: helium_env 1191 REAL(KIND=dp), INTENT(IN) :: initkT, solute_radius 1192 1193 CHARACTER(len=*), PARAMETER :: routineN = 'helium_coord_init', & 1194 routineP = moduleN//':'//routineN 1195 REAL(KIND=dp), PARAMETER :: minHeHedst = 5.669177966_dp 1196 1197 INTEGER :: ia, ib, ic, id, iter, k 1198 LOGICAL :: invalidpos 1199 REAL(KIND=dp) :: minHeHedsttmp, r1, r2 1200 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: centroids 1201 REAL(KIND=dp), DIMENSION(3) :: cvek, rvek, tvek 1202 1203 !corresponds to three angstrom (roughly first maximum of He-He-rdf) 1204 minHeHedsttmp = minHeHedst 1205 1206 DO k = 1, SIZE(helium_env) 1207 IF (helium_env(k)%helium%solute_present) THEN 1208 cvek(:) = helium_env(k)%helium%center(:) 1209 CALL helium_pbc(helium_env(k)%helium, cvek) 1210 END IF 1211 1212 ALLOCATE (centroids(3, helium_env(k)%helium%atoms)) 1213 IF (helium_env(k)%helium%periodic) THEN 1214 DO ia = 1, helium_env(k)%helium%atoms 1215 invalidpos = .TRUE. 1216 iter = 0 1217 DO WHILE (invalidpos) 1218 iter = iter + 1 1219 invalidpos = .FALSE. 1220 ! if sampling fails to often, reduce he he criterion 1221 !CS TODO: 1222 !minHeHedsttmp = 0.90_dp**(iter/100)*minHeHedst 1223 minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst 1224 DO ic = 1, 3 1225 r1 = next_random_number(helium_env(k)%helium%rng_stream_uniform) 1226 r1 = 2.0_dp*r1 - 1.0_dp 1227 r1 = r1*helium_env(k)%helium%cell_size 1228 centroids(ic, ia) = r1 1229 END DO 1230 ! check if helium is outside of cell 1231 tvek(:) = centroids(:, ia) 1232 CALL helium_pbc(helium_env(k)%helium, tvek(:)) 1233 rvek(:) = tvek(:) - centroids(:, ia) 1234 r2 = DOT_PRODUCT(rvek, rvek) 1235 IF (r2 > 1.0_dp*10.0_dp**(-6)) THEN 1236 invalidpos = .TRUE. 1237 ELSE 1238 ! check for helium-helium collision 1239 DO id = 1, ia - 1 1240 rvek = centroids(:, ia) - centroids(:, id) 1241 CALL helium_pbc(helium_env(k)%helium, rvek) 1242 r2 = DOT_PRODUCT(rvek, rvek) 1243 IF (r2 < minHeHedsttmp**2) THEN 1244 invalidpos = .TRUE. 1245 EXIT 1246 END IF 1247 END DO 1248 END IF 1249 IF (.NOT. invalidpos) THEN 1250 ! check if centroid collides with molecule 1251 IF (helium_env(k)%helium%solute_present) THEN 1252 rvek(:) = (cvek(:) - centroids(:, ia)) 1253 r2 = DOT_PRODUCT(rvek, rvek) 1254 IF (r2 <= solute_radius**2) invalidpos = .TRUE. 1255 END IF 1256 END IF 1257 END DO 1258 END DO 1259 ! do thermal gaussian delocalization of hot start 1260 IF (initkT > 0.0_dp) THEN 1261 CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkT) 1262 ELSE 1263 DO ia = 1, helium_env(k)%helium%atoms 1264 DO ib = 1, helium_env(k)%helium%beads 1265 helium_env(k)%helium%pos(:, ia, ib) = centroids(:, ia) 1266 END DO 1267 END DO 1268 END IF 1269 ! apply PBC to bead coords 1270 DO ia = 1, helium_env(k)%helium%atoms 1271 DO ib = 1, helium_env(k)%helium%beads 1272 CALL helium_pbc(helium_env(k)%helium, helium_env(k)%helium%pos(:, ia, ib)) 1273 ! check if bead collides with molecule 1274 IF (helium_env(k)%helium%solute_present) THEN 1275 rvek(:) = (cvek(:) - helium_env(k)%helium%pos(:, ia, ib)) 1276 r2 = DOT_PRODUCT(rvek, rvek) 1277 IF (r2 <= solute_radius**2) THEN 1278 r1 = SQRT(r2) 1279 helium_env(k)%helium%pos(:, ia, ib) = & 1280 cvek(:) + solute_radius/r1*rvek(:) 1281 END IF 1282 END IF 1283 END DO 1284 END DO 1285 ELSE 1286 DO ia = 1, helium_env(k)%helium%atoms 1287 iter = 0 1288 invalidpos = .TRUE. 1289 DO WHILE (invalidpos) 1290 invalidpos = .FALSE. 1291 iter = iter + 1 1292 ! if sampling fails to often, reduce he he criterion 1293 minHeHedsttmp = 0.90_dp**MIN(0, iter - 2)*minHeHedst 1294 DO ic = 1, 3 1295 rvek(ic) = next_random_number(helium_env(k)%helium%rng_stream_uniform) 1296 rvek(ic) = 2.0_dp*rvek(ic) - 1.0_dp 1297 rvek(ic) = rvek(ic)*helium_env(k)%helium%droplet_radius 1298 END DO 1299 centroids(:, ia) = rvek(:) 1300 ! check if helium is outside of the droplet 1301 r2 = DOT_PRODUCT(rvek, rvek) 1302 IF (r2 > helium_env(k)%helium%droplet_radius**2) THEN 1303 invalidpos = .TRUE. 1304 ELSE 1305 ! check for helium-helium collision 1306 DO id = 1, ia - 1 1307 rvek = centroids(:, ia) - centroids(:, id) 1308 r2 = DOT_PRODUCT(rvek, rvek) 1309 IF (r2 < minHeHedsttmp**2) THEN 1310 invalidpos = .TRUE. 1311 EXIT 1312 END IF 1313 END DO 1314 END IF 1315 IF (.NOT. invalidpos) THEN 1316 ! make sure the helium does not collide with the solute 1317 IF (helium_env(k)%helium%solute_present) THEN 1318 rvek(:) = (cvek(:) - centroids(:, ia)) 1319 r2 = DOT_PRODUCT(rvek, rvek) 1320 IF (r2 <= solute_radius**2) invalidpos = .TRUE. 1321 END IF 1322 END IF 1323 END DO 1324 END DO 1325 ! do thermal gaussian delocalization of hot start 1326 IF (initkT > 0.0_dp) THEN 1327 CALL helium_thermal_gaussian_beads_init(helium_env(k)%helium, centroids, initkT) 1328 ELSE 1329 DO ia = 1, helium_env(k)%helium%atoms 1330 DO ib = 1, helium_env(k)%helium%beads 1331 helium_env(k)%helium%pos(:, ia, ib) = centroids(:, ia) 1332 END DO 1333 END DO 1334 END IF 1335 DO ia = 1, helium_env(k)%helium%atoms 1336 DO ib = 1, helium_env(k)%helium%beads 1337 ! Make sure, that nothing lies outside the droplet radius 1338 r1 = DOT_PRODUCT(helium_env(k)%helium%pos(:, ia, ib), & 1339 helium_env(k)%helium%pos(:, ia, ib)) 1340 IF (r1 > helium_env(k)%helium%droplet_radius**2) THEN 1341 r1 = SQRT(r1) 1342 helium_env(k)%helium%pos(:, ia, ib) = & 1343 helium_env(k)%helium%droplet_radius/r1* & 1344 helium_env(k)%helium%pos(:, ia, ib) 1345 ELSE IF (helium_env(k)%helium%solute_present) THEN 1346 IF (r1 < solute_radius**2) THEN 1347 !make sure that nothing lies within the molecule 1348 r1 = SQRT(r1) 1349 helium_env(k)%helium%pos(:, ia, ib) = & 1350 solute_radius/r1* & 1351 helium_env(k)%helium%pos(:, ia, ib) 1352 END IF 1353 END IF 1354 ! transfer to position around actual center of droplet 1355 helium_env(k)%helium%pos(:, ia, ib) = & 1356 helium_env(k)%helium%pos(:, ia, ib) + & 1357 helium_env(k)%helium%center(:) 1358 END DO 1359 END DO 1360 END IF 1361 helium_env(k)%helium%work = helium_env(k)%helium%pos 1362 DEALLOCATE (centroids) 1363 END DO 1364 1365 RETURN 1366 END SUBROUTINE helium_coord_init 1367 1368! ************************************************************************************************** 1369!> \brief ... 1370!> \param helium_env ... 1371!> \param centroids ... 1372!> \param kbT ... 1373! ************************************************************************************************** 1374 SUBROUTINE helium_thermal_gaussian_beads_init(helium_env, centroids, kbT) 1375 1376 TYPE(helium_solvent_type), POINTER :: helium_env 1377 REAL(KIND=dp), DIMENSION(3, helium_env%atoms), & 1378 INTENT(IN) :: centroids 1379 REAL(KIND=dp), INTENT(IN) :: kbT 1380 1381 INTEGER :: i, iatom, idim, imode, j, p 1382 REAL(KIND=dp) :: invsqrtp, omega, pip, rand, sqrt2p, & 1383 sqrtp, twopip, variance 1384 REAL(KIND=dp), & 1385 DIMENSION(helium_env%beads, helium_env%beads) :: u2x 1386 REAL(KIND=dp), DIMENSION(helium_env%beads) :: nmhecoords 1387 1388 p = helium_env%beads 1389 1390 sqrt2p = SQRT(2.0_dp/REAL(p, dp)) 1391 twopip = twopi/REAL(p, dp) 1392 pip = pi/REAL(p, dp) 1393 sqrtp = SQRT(REAL(p, dp)) 1394 invsqrtp = 1.0_dp/SQRT(REAL(p, dp)) 1395 1396 ! set up normal mode backtransform matrix 1397 u2x(:, :) = 0.0_dp 1398 u2x(:, 1) = invsqrtp 1399 DO i = 2, p/2 + 1 1400 DO j = 1, p 1401 u2x(j, i) = sqrt2p*COS(twopip*(i - 1)*(j - 1)) 1402 END DO 1403 END DO 1404 DO i = p/2 + 2, p 1405 DO j = 1, p 1406 u2x(j, i) = sqrt2p*SIN(twopip*(i - 1)*(j - 1)) 1407 END DO 1408 END DO 1409 IF (MOD(p, 2) == 0) THEN 1410 DO i = 1, p - 1, 2 1411 u2x(i, p/2 + 1) = invsqrtp 1412 u2x(i + 1, p/2 + 1) = -1.0_dp*invsqrtp 1413 END DO 1414 END IF 1415 1416 DO iatom = 1, helium_env%atoms 1417 DO idim = 1, 3 1418 nmhecoords(1) = sqrtp*centroids(idim, iatom) 1419 DO imode = 2, p 1420 omega = 2.0_dp*p*kbT*SIN((imode - 1)*pip) 1421 variance = kbT*p/(helium_env%he_mass_au*omega**2) 1422 rand = next_random_number(helium_env%rng_stream_gaussian) 1423 nmhecoords(imode) = rand*SQRT(variance) 1424 END DO 1425 helium_env%pos(idim, iatom, 1:p) = MATMUL(u2x, nmhecoords) 1426 END DO 1427 END DO 1428 1429 END SUBROUTINE helium_thermal_gaussian_beads_init 1430 1431! *************************************************************************** 1432!> \brief Restore coordinates from the input structure. 1433!> \param helium_env ... 1434!> \date 2009-11-09 1435!> \par History 1436!> 2010-07-22 accomodate additional cpus in the runtime wrt the 1437!> restart [lwalewski] 1438!> 2016-07-14 Modified to work with independent helium_env 1439!> [cschran] 1440!> \author Lukasz Walewski 1441! ************************************************************************************************** 1442 SUBROUTINE helium_coord_restore(helium_env) 1443 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1444 1445 CHARACTER(len=*), PARAMETER :: routineN = 'helium_coord_restore', & 1446 routineP = moduleN//':'//routineN 1447 1448 CHARACTER(len=default_string_length) :: err_str, stmp 1449 INTEGER :: actlen, i, k, msglen, num_env_restart, & 1450 off, offset 1451 LOGICAL, DIMENSION(:, :, :), POINTER :: m 1452 REAL(KIND=dp), DIMENSION(:), POINTER :: message 1453 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: f 1454 TYPE(cp_logger_type), POINTER :: logger 1455 1456 NULLIFY (logger) 1457 logger => cp_get_default_logger() 1458 1459 ! assign the pointer to the memory location of the input structure, where 1460 ! the coordinates are stored 1461 NULLIFY (message) 1462 CALL section_vals_val_get(helium_env(1)%helium%input, & 1463 "MOTION%PINT%HELIUM%COORD%_DEFAULT_KEYWORD_", & 1464 r_vals=message) 1465 1466 ! check that the number of values in the input match the current runtime 1467 actlen = SIZE(message) 1468 num_env_restart = actlen/helium_env(1)%helium%atoms/helium_env(1)%helium%beads/3 1469 1470 IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN 1471 err_str = "Reading bead coordinates from the input file." 1472 CALL helium_write_line(err_str) 1473 err_str = "Number of environments in the restart...: '" 1474 stmp = "" 1475 WRITE (stmp, *) num_env_restart 1476 err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'." 1477 CALL helium_write_line(err_str) 1478 err_str = "Number of current run time environments.: '" 1479 stmp = "" 1480 WRITE (stmp, *) helium_env(1)%helium%num_env 1481 err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'." 1482 CALL helium_write_line(err_str) 1483 err_str = "Missmatch between number of bead coord. in input file and helium environments." 1484 CPABORT(err_str) 1485 ELSE 1486 CALL helium_write_line("Bead coordinates read from the input file.") 1487 1488 offset = 0 1489 DO i = 1, logger%para_env%mepos 1490 offset = offset + helium_env(1)%env_all(i) 1491 END DO 1492 1493 ! distribute coordinates over processors (no message passing) 1494 DO k = 1, SIZE(helium_env) 1495 msglen = helium_env(k)%helium%atoms*helium_env(k)%helium%beads*3 1496 off = msglen*MOD(offset + k - 1, num_env_restart) 1497 NULLIFY (m, f) 1498 ALLOCATE (m(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads)) 1499 ALLOCATE (f(3, helium_env(k)%helium%atoms, helium_env(k)%helium%beads)) 1500 m(:, :, :) = .TRUE. 1501 f(:, :, :) = 0.0_dp 1502 helium_env(k)%helium%pos(:, :, 1:helium_env(k)%helium%beads) = UNPACK(message(off + 1:off + msglen), MASK=m, FIELD=f) 1503 DEALLOCATE (f, m) 1504 END DO 1505 1506 END IF 1507 1508 NULLIFY (message) 1509 1510 RETURN 1511 END SUBROUTINE helium_coord_restore 1512 1513! *************************************************************************** 1514!> \brief Initialize forces exerted on the solute 1515!> \param helium_env ... 1516!> \date 2009-11-10 1517!> \par History 1518!> 2016-07-14 Modified to work with independent helium_env [cschran] 1519!> \author Lukasz Walewski 1520! ************************************************************************************************** 1521 SUBROUTINE helium_force_init(helium_env) 1522 1523 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1524 1525 CHARACTER(len=*), PARAMETER :: routineN = 'helium_force_init', & 1526 routineP = moduleN//':'//routineN 1527 1528 INTEGER :: k 1529 1530 DO k = 1, SIZE(helium_env) 1531 IF (helium_env(k)%helium%solute_present) THEN 1532 helium_env(k)%helium%force_avrg(:, :) = 0.0_dp 1533 helium_env(k)%helium%force_inst(:, :) = 0.0_dp 1534 END IF 1535 END DO 1536 1537 RETURN 1538 END SUBROUTINE helium_force_init 1539 1540! *************************************************************************** 1541!> \brief Restore forces from the input structure to the runtime environment. 1542!> \param helium_env ... 1543!> \date 2009-11-10 1544!> \par History 1545!> 2016-07-14 Modified to work with independent helium_env [cschran] 1546!> \author Lukasz Walewski 1547! ************************************************************************************************** 1548 SUBROUTINE helium_force_restore(helium_env) 1549 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1550 1551 CHARACTER(len=*), PARAMETER :: routineN = 'helium_force_restore', & 1552 routineP = moduleN//':'//routineN 1553 1554 CHARACTER(len=default_string_length) :: err_str, stmp 1555 INTEGER :: actlen, k, msglen 1556 LOGICAL, DIMENSION(:, :), POINTER :: m 1557 REAL(KIND=dp), DIMENSION(:), POINTER :: message 1558 REAL(KIND=dp), DIMENSION(:, :), POINTER :: f 1559 1560! assign the pointer to the memory location of the input structure, where 1561! the forces are stored 1562 1563 NULLIFY (message) 1564 CALL section_vals_val_get(helium_env(1)%helium%input, & 1565 "MOTION%PINT%HELIUM%FORCE%_DEFAULT_KEYWORD_", & 1566 r_vals=message) 1567 1568 ! check if the destination array has correct size 1569 msglen = helium_env(1)%helium%solute_atoms*helium_env(1)%helium%solute_beads*3 1570 actlen = SIZE(helium_env(1)%helium%force_avrg) 1571 err_str = "Invalid size of helium%force_avrg array: actual '" 1572 stmp = "" 1573 WRITE (stmp, *) actlen 1574 err_str = TRIM(ADJUSTL(err_str))// & 1575 TRIM(ADJUSTL(stmp))//"' but expected '" 1576 stmp = "" 1577 WRITE (stmp, *) msglen 1578 IF (actlen /= msglen) THEN 1579 err_str = TRIM(ADJUSTL(err_str))// & 1580 TRIM(ADJUSTL(stmp))//"'." 1581 CPABORT(err_str) 1582 END IF 1583 1584 ! restore forces on all processors (no message passing) 1585 NULLIFY (m, f) 1586 ALLOCATE (m(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3)) 1587 ALLOCATE (f(helium_env(1)%helium%solute_beads, helium_env(1)%helium%solute_atoms*3)) 1588 m(:, :) = .TRUE. 1589 f(:, :) = 0.0_dp 1590 DO k = 1, SIZE(helium_env) 1591 helium_env(k)%helium%force_avrg(:, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f) 1592 helium_env(k)%helium%force_inst(:, :) = 0.0_dp 1593 END DO 1594 DEALLOCATE (f, m) 1595 1596 CALL helium_write_line("Forces on the solute read from the input file.") 1597 1598 NULLIFY (message) 1599 1600 RETURN 1601 END SUBROUTINE helium_force_restore 1602 1603! *************************************************************************** 1604!> \brief Initialize the permutation state. 1605!> \param helium_env ... 1606!> \date 2009-11-05 1607!> \par History 1608!> 2016-07-14 Modified to work with independent helium_env [cschran] 1609!> \author Lukasz Walewski 1610!> \note Assign the identity permutation at each processor. Inverse 1611!> permutation array gets assigned as well. 1612! ************************************************************************************************** 1613 SUBROUTINE helium_perm_init(helium_env) 1614 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1615 1616 CHARACTER(len=*), PARAMETER :: routineN = 'helium_perm_init', & 1617 routineP = moduleN//':'//routineN 1618 1619 INTEGER :: ia, k 1620 1621 DO k = 1, SIZE(helium_env) 1622 DO ia = 1, helium_env(k)%helium%atoms 1623 helium_env(k)%helium%permutation(ia) = ia 1624 helium_env(k)%helium%iperm(ia) = ia 1625 END DO 1626 END DO 1627 1628 RETURN 1629 END SUBROUTINE helium_perm_init 1630 1631! *************************************************************************** 1632!> \brief Restore permutation state from the input structre. 1633!> \param helium_env ... 1634!> \date 2009-11-05 1635!> \par History 1636!> 2010-07-22 accomodate additional cpus in the runtime wrt the 1637!> restart [lwalewski] 1638!> 2016-07-14 Modified to work with independent helium_env [cschran] 1639!> \author Lukasz Walewski 1640!> \note Transfer permutation state from the input tree to the runtime 1641!> data structures on each processor. Inverse permutation array is 1642!> recalculated according to the restored permutation state. 1643! ************************************************************************************************** 1644 SUBROUTINE helium_perm_restore(helium_env) 1645 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1646 1647 CHARACTER(len=*), PARAMETER :: routineN = 'helium_perm_restore', & 1648 routineP = moduleN//':'//routineN 1649 1650 CHARACTER(len=default_string_length) :: err_str, stmp 1651 INTEGER :: actlen, i, ia, ic, k, msglen, & 1652 num_env_restart, off, offset 1653 INTEGER, DIMENSION(:), POINTER :: message 1654 TYPE(cp_logger_type), POINTER :: logger 1655 1656 NULLIFY (logger) 1657 logger => cp_get_default_logger() 1658 1659 ! assign the pointer to the memory location of the input structure, where 1660 ! the permutation state is stored 1661 NULLIFY (message) 1662 CALL section_vals_val_get(helium_env(1)%helium%input, & 1663 "MOTION%PINT%HELIUM%PERM%_DEFAULT_KEYWORD_", & 1664 i_vals=message) 1665 1666 ! check the number of environments presumably stored in the restart 1667 actlen = SIZE(message) 1668 num_env_restart = actlen/helium_env(1)%helium%atoms 1669 1670!TODO maybe add some sanity checks here: 1671! is num_env_restart integer ? 1672 1673 IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN 1674 err_str = "Reading permutation state from the input file." 1675 CALL helium_write_line(err_str) 1676 err_str = "Number of environments in the restart...: '" 1677 stmp = "" 1678 WRITE (stmp, *) num_env_restart 1679 err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'." 1680 CALL helium_write_line(err_str) 1681 err_str = "Number of current run time environments.: '" 1682 stmp = "" 1683 WRITE (stmp, *) helium_env(1)%helium%num_env 1684 err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'." 1685 CALL helium_write_line(err_str) 1686 err_str = "Missmatch between number of perm. states in input file and helium environments." 1687 CPABORT(err_str) 1688 ELSE 1689 CALL helium_write_line("Permutation state read from the input file.") 1690 1691 ! distribute permutation state over processors 1692 offset = 0 1693 DO i = 1, logger%para_env%mepos 1694 offset = offset + helium_env(1)%env_all(i) 1695 END DO 1696 1697 DO k = 1, SIZE(helium_env) 1698 msglen = helium_env(k)%helium%atoms 1699 off = msglen*MOD(k - 1 + offset, num_env_restart) 1700 helium_env(k)%helium%permutation(:) = message(off + 1:off + msglen) 1701 END DO 1702 END IF 1703 1704 ! recalculate the inverse permutation array 1705 DO k = 1, SIZE(helium_env) 1706 helium_env(k)%helium%iperm(:) = 0 1707 ic = 0 1708 DO ia = 1, msglen 1709 IF ((helium_env(k)%helium%permutation(ia) > 0) .AND. (helium_env(k)%helium%permutation(ia) <= msglen)) THEN 1710 helium_env(k)%helium%iperm(helium_env(k)%helium%permutation(ia)) = ia 1711 ic = ic + 1 1712 END IF 1713 END DO 1714 err_str = "Invalid HELIUM%PERM state: some numbers not within (1," 1715 stmp = "" 1716 WRITE (stmp, *) msglen 1717 IF (ic /= msglen) THEN 1718 err_str = TRIM(ADJUSTL(err_str))// & 1719 TRIM(ADJUSTL(stmp))//")." 1720 CPABORT(err_str) 1721 END IF 1722 END DO 1723 NULLIFY (message) 1724 1725 RETURN 1726 END SUBROUTINE helium_perm_restore 1727 1728! *************************************************************************** 1729!> \brief Restore averages from the input structure 1730!> \param helium_env ... 1731!> \date 2014-06-25 1732!> \par History 1733!> 2016-07-14 Modified to work with independent helium_env [cschran] 1734!> \author Lukasz Walewski 1735! ************************************************************************************************** 1736 SUBROUTINE helium_averages_restore(helium_env) 1737 1738 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1739 1740 CHARACTER(len=*), PARAMETER :: routineN = 'helium_averages_restore', & 1741 routineP = moduleN//':'//routineN 1742 1743 INTEGER :: i, k, msglen, num_env_restart, off, & 1744 offset 1745 LOGICAL :: explicit 1746 REAL(KIND=dp), DIMENSION(:), POINTER :: message 1747 TYPE(cp_logger_type), POINTER :: logger 1748 1749 NULLIFY (logger) 1750 logger => cp_get_default_logger() 1751 1752 offset = 0 1753 DO i = 1, logger%para_env%mepos 1754 offset = offset + helium_env(1)%env_all(i) 1755 END DO 1756 1757 ! restore projected area 1758 CALL section_vals_val_get(helium_env(1)%helium%input, & 1759 "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", & 1760 explicit=explicit) 1761 IF (explicit) THEN 1762 NULLIFY (message) 1763 CALL section_vals_val_get(helium_env(1)%helium%input, & 1764 "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA", & 1765 r_vals=message) 1766 num_env_restart = SIZE(message)/3 ! apparent number of environments 1767 msglen = 3 1768 DO k = 1, SIZE(helium_env) 1769 off = msglen*MOD(offset + k - 1, num_env_restart) 1770 helium_env(k)%helium%proarea%rstr(:) = message(off + 1:off + msglen) 1771 END DO 1772 ELSE 1773 DO k = 1, SIZE(helium_env) 1774 helium_env(k)%helium%proarea%rstr(:) = 0.0_dp 1775 END DO 1776 END IF 1777 1778 ! restore projected area squared 1779 CALL section_vals_val_get(helium_env(1)%helium%input, & 1780 "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", & 1781 explicit=explicit) 1782 IF (explicit) THEN 1783 NULLIFY (message) 1784 CALL section_vals_val_get(helium_env(1)%helium%input, & 1785 "MOTION%PINT%HELIUM%AVERAGES%PROJECTED_AREA_2", & 1786 r_vals=message) 1787 num_env_restart = SIZE(message)/3 ! apparent number of environments 1788 msglen = 3 1789 DO k = 1, SIZE(helium_env) 1790 off = msglen*MOD(offset + k - 1, num_env_restart) 1791 helium_env(k)%helium%prarea2%rstr(:) = message(off + 1:off + msglen) 1792 END DO 1793 ELSE 1794 DO k = 1, SIZE(helium_env) 1795 helium_env(k)%helium%prarea2%rstr(:) = 0.0_dp 1796 END DO 1797 END IF 1798 1799 ! restore winding number squared 1800 CALL section_vals_val_get(helium_env(1)%helium%input, & 1801 "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", & 1802 explicit=explicit) 1803 IF (explicit) THEN 1804 NULLIFY (message) 1805 CALL section_vals_val_get(helium_env(1)%helium%input, & 1806 "MOTION%PINT%HELIUM%AVERAGES%WINDING_NUMBER_2", & 1807 r_vals=message) 1808 num_env_restart = SIZE(message)/3 ! apparent number of environments 1809 msglen = 3 1810 DO k = 1, SIZE(helium_env) 1811 off = msglen*MOD(offset + k - 1, num_env_restart) 1812 helium_env(k)%helium%wnmber2%rstr(:) = message(off + 1:off + msglen) 1813 END DO 1814 ELSE 1815 DO k = 1, SIZE(helium_env) 1816 helium_env(k)%helium%wnmber2%rstr(:) = 0.0_dp 1817 END DO 1818 END IF 1819 1820 ! restore moment of inertia 1821 CALL section_vals_val_get(helium_env(1)%helium%input, & 1822 "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", & 1823 explicit=explicit) 1824 IF (explicit) THEN 1825 NULLIFY (message) 1826 CALL section_vals_val_get(helium_env(1)%helium%input, & 1827 "MOTION%PINT%HELIUM%AVERAGES%MOMENT_OF_INERTIA", & 1828 r_vals=message) 1829 num_env_restart = SIZE(message)/3 ! apparent number of environments 1830 msglen = 3 1831 DO k = 1, SIZE(helium_env) 1832 off = msglen*MOD(offset + k - 1, num_env_restart) 1833 helium_env(k)%helium%mominer%rstr(:) = message(off + 1:off + msglen) 1834 END DO 1835 ELSE 1836 DO k = 1, SIZE(helium_env) 1837 helium_env(k)%helium%mominer%rstr(:) = 0.0_dp 1838 END DO 1839 END IF 1840 1841 IF (helium_env(1)%helium%rdf_present) THEN 1842 CALL helium_rdf_restore(helium_env) 1843 END IF 1844 1845 IF (helium_env(1)%helium%rho_present) THEN 1846 ! restore densities 1847 CALL helium_rho_restore(helium_env) 1848 END IF 1849 1850 ! get the weighting factor 1851 DO k = 1, SIZE(helium_env) 1852 CALL section_vals_val_get( & 1853 helium_env(k)%helium%input, & 1854 "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", & 1855 i_val=helium_env(k)%helium%averages_iweight) 1856 1857 ! set the flag indicating whether the averages have been restarted 1858 CALL section_vals_val_get( & 1859 helium_env(k)%helium%input, & 1860 "EXT_RESTART%RESTART_HELIUM_AVERAGES", & 1861 l_val=helium_env(k)%helium%averages_restarted) 1862 END DO 1863 1864 RETURN 1865 END SUBROUTINE helium_averages_restore 1866 1867! *************************************************************************** 1868!> \brief Create RNG streams and initialize their state. 1869!> \param helium_env ... 1870!> \date 2009-11-04 1871!> \par History 1872!> 2016-07-14 Modified to work with independent helium_env [cschran] 1873!> \author Lukasz Walewski 1874!> \note TODO: This function shouldn't create (allocate) objects! Only 1875!> initialization, i.e. setting the seed values etc, should be done 1876!> here, allocation should be moved to helium_create 1877! ************************************************************************************************** 1878 SUBROUTINE helium_rng_init(helium_env) 1879 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1880 1881 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rng_init', & 1882 routineP = moduleN//':'//routineN 1883 1884 INTEGER :: helium_seed, i, offset 1885 REAL(KIND=dp), DIMENSION(3, 2) :: initial_seed 1886 TYPE(cp_logger_type), POINTER :: logger 1887 TYPE(rng_stream_p_type), DIMENSION(:), POINTER :: gaussian_array, uniform_array 1888 TYPE(rng_stream_type), POINTER :: next_rngs, prev_rngs 1889 1890 NULLIFY (logger) 1891 logger => cp_get_default_logger() 1892 IF (logger%para_env%ionode) THEN 1893 CALL section_vals_val_get(helium_env(1)%helium%input, & 1894 "MOTION%PINT%HELIUM%RNG_SEED", & 1895 i_val=helium_seed) 1896 END IF 1897 CALL mp_bcast(helium_seed, & 1898 logger%para_env%source, & 1899 helium_env(1)%comm) 1900 initial_seed(:, :) = REAL(helium_seed, dp) 1901 1902 ALLOCATE (uniform_array(helium_env(1)%helium%num_env), & 1903 gaussian_array(helium_env(1)%helium%num_env)) 1904 DO i = 1, helium_env(1)%helium%num_env 1905 NULLIFY (uniform_array(i)%stream, gaussian_array(i)%stream) 1906 END DO 1907 NULLIFY (prev_rngs, next_rngs) 1908 1909 ! Create num_env RNG streams on processor all processors 1910 ! and distribute them so, that each processor gets unique 1911 ! RN sequences for his helium environments 1912 ! COMMENT: rng_stream can not be used with mp_bcast 1913 1914 CALL create_rng_stream(prev_rngs, & 1915 name="helium_rns_uniform", & 1916 distribution_type=UNIFORM, & 1917 extended_precision=.TRUE., & 1918 seed=initial_seed) 1919 uniform_array(1)%stream => prev_rngs 1920 1921 CALL create_rng_stream(next_rngs, & 1922 name="helium_rns_gaussian", & 1923 last_rng_stream=prev_rngs, & 1924 distribution_type=GAUSSIAN, & 1925 extended_precision=.TRUE.) 1926 gaussian_array(1)%stream => next_rngs 1927 NULLIFY (prev_rngs) 1928 prev_rngs => next_rngs 1929 NULLIFY (next_rngs) 1930 1931 DO i = 2, helium_env(1)%helium%num_env 1932 CALL create_rng_stream(next_rngs, & 1933 name="helium_rns_uniform", & 1934 last_rng_stream=prev_rngs, & 1935 distribution_type=UNIFORM, & 1936 extended_precision=.TRUE.) 1937 uniform_array(i)%stream => next_rngs 1938 prev_rngs => next_rngs 1939 NULLIFY (next_rngs) 1940 1941 CALL create_rng_stream(next_rngs, & 1942 name="helium_rns_gaussian", & 1943 last_rng_stream=prev_rngs, & 1944 distribution_type=GAUSSIAN, & 1945 extended_precision=.TRUE.) 1946 gaussian_array(i)%stream => next_rngs 1947 prev_rngs => next_rngs 1948 NULLIFY (next_rngs) 1949 1950 END DO 1951 1952 NULLIFY (prev_rngs) 1953 1954 offset = 0 1955 DO i = 1, logger%para_env%mepos 1956 offset = offset + helium_env(1)%env_all(i) 1957 END DO 1958 1959 IF (ASSOCIATED(helium_env)) THEN 1960 DO i = 1, SIZE(helium_env) 1961 NULLIFY (helium_env(i)%helium%rng_stream_uniform, & 1962 helium_env(i)%helium%rng_stream_gaussian) 1963 helium_env(i)%helium%rng_stream_uniform => uniform_array(offset + i)%stream 1964 helium_env(i)%helium%rng_stream_gaussian => gaussian_array(offset + i)%stream 1965 END DO 1966 END IF 1967 1968 DO i = 1, helium_env(1)%helium%num_env 1969 IF (i .LE. offset .OR. i .GT. offset + SIZE(helium_env)) THEN 1970 CALL delete_rng_stream(uniform_array(i)%stream) 1971 CALL delete_rng_stream(gaussian_array(i)%stream) 1972 END IF 1973 NULLIFY (uniform_array(i)%stream) 1974 NULLIFY (gaussian_array(i)%stream) 1975 END DO 1976 1977 DEALLOCATE (uniform_array) 1978 DEALLOCATE (gaussian_array) 1979 1980 RETURN 1981 END SUBROUTINE helium_rng_init 1982 1983! *************************************************************************** 1984!> \brief Restore RNG state from the input structure. 1985!> \param helium_env ... 1986!> \date 2009-11-04 1987!> \par History 1988!> 2010-07-22 Create new rng streams if more cpus available in the 1989!> runtime than in the restart [lwalewski] 1990!> 2016-04-18 Modified for independet number of helium_env [cschran] 1991!> \author Lukasz Walewski 1992! ************************************************************************************************** 1993 SUBROUTINE helium_rng_restore(helium_env) 1994 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 1995 1996 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rng_restore', & 1997 routineP = moduleN//':'//routineN 1998 1999 CHARACTER(len=default_string_length) :: err_str, stmp 2000 INTEGER :: actlen, i, k, msglen, num_env_restart, & 2001 off, offset 2002 LOGICAL :: lbf 2003 LOGICAL, DIMENSION(3, 2) :: m 2004 REAL(KIND=dp) :: bf, bu 2005 REAL(KIND=dp), DIMENSION(3, 2) :: bg, cg, f, ig 2006 REAL(KIND=dp), DIMENSION(:), POINTER :: message 2007 TYPE(cp_logger_type), POINTER :: logger 2008 2009 NULLIFY (logger) 2010 logger => cp_get_default_logger() 2011 2012 ! assign the pointer to the memory location of the input structure 2013 ! where the RNG state is stored 2014 NULLIFY (message) 2015 CALL section_vals_val_get(helium_env(1)%helium%input, & 2016 "MOTION%PINT%HELIUM%RNG_STATE%_DEFAULT_KEYWORD_", & 2017 r_vals=message) 2018 2019 ! check the number of environments presumably stored in the restart 2020 actlen = SIZE(message) 2021 num_env_restart = actlen/40 2022 2023 ! check, if RNG restart has the same dimension as helium%num_env 2024 IF (num_env_restart .NE. helium_env(1)%helium%num_env) THEN 2025 err_str = "Reading RNG state from the input file." 2026 CALL helium_write_line(err_str) 2027 err_str = "Number of environments in the restart...: '" 2028 stmp = "" 2029 WRITE (stmp, *) num_env_restart 2030 err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'." 2031 CALL helium_write_line(err_str) 2032 err_str = "Number of current run time environments.: '" 2033 stmp = "" 2034 WRITE (stmp, *) helium_env(1)%helium%num_env 2035 err_str = TRIM(ADJUSTL(err_str))//TRIM(ADJUSTL(stmp))//"'." 2036 CALL helium_write_line(err_str) 2037 err_str = "Missmatch between number of RNG states in input file and helium environments." 2038 CPABORT(err_str) 2039 ELSE 2040 CALL helium_write_line("RNG state read from the input file.") 2041 2042 ! unpack the buffer at each processor, set RNG state 2043 offset = 0 2044 DO i = 1, logger%para_env%mepos 2045 offset = offset + helium_env(1)%env_all(i) 2046 END DO 2047 2048 DO k = 1, SIZE(helium_env) 2049 msglen = 40 2050 off = msglen*(offset + k - 1) 2051 m(:, :) = .TRUE. 2052 f(:, :) = 0.0_dp 2053 bg(:, :) = UNPACK(message(off + 1:off + 6), MASK=m, FIELD=f) 2054 cg(:, :) = UNPACK(message(off + 7:off + 12), MASK=m, FIELD=f) 2055 ig(:, :) = UNPACK(message(off + 13:off + 18), MASK=m, FIELD=f) 2056 bf = message(off + 19) 2057 bu = message(off + 20) 2058 IF (bf .GT. 0) THEN 2059 lbf = .TRUE. 2060 ELSE 2061 lbf = .FALSE. 2062 END IF 2063 CALL set_rng_stream(helium_env(k)%helium%rng_stream_uniform, bg=bg, cg=cg, ig=ig, & 2064 buffer=bu, buffer_filled=lbf) 2065 bg(:, :) = UNPACK(message(off + 21:off + 26), MASK=m, FIELD=f) 2066 cg(:, :) = UNPACK(message(off + 27:off + 32), MASK=m, FIELD=f) 2067 ig(:, :) = UNPACK(message(off + 33:off + 38), MASK=m, FIELD=f) 2068 bf = message(off + 39) 2069 bu = message(off + 40) 2070 IF (bf .GT. 0) THEN 2071 lbf = .TRUE. 2072 ELSE 2073 lbf = .FALSE. 2074 END IF 2075 CALL set_rng_stream(helium_env(k)%helium%rng_stream_gaussian, bg=bg, cg=cg, ig=ig, & 2076 buffer=bu, buffer_filled=lbf) 2077 END DO 2078 END IF 2079 2080 NULLIFY (message) 2081 2082 RETURN 2083 END SUBROUTINE helium_rng_restore 2084 2085! *************************************************************************** 2086!> \brief Create the RDF related data structures and initialize 2087!> \param helium ... 2088!> \date 2014-02-25 2089!> \par History 2090!> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran] 2091!> \author Lukasz Walewski 2092! ************************************************************************************************** 2093 SUBROUTINE helium_rdf_init(helium) 2094 2095 TYPE(helium_solvent_type), POINTER :: helium 2096 2097 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rdf_init', & 2098 routineP = moduleN//':'//routineN 2099 2100 CHARACTER(len=2*default_string_length) :: err_str, stmp 2101 INTEGER :: ii, ij 2102 LOGICAL :: explicit 2103 TYPE(cp_logger_type), POINTER :: logger 2104 2105 ! read parameters 2106 NULLIFY (logger) 2107 logger => cp_get_default_logger() 2108 CALL section_vals_val_get( & 2109 helium%input, & 2110 "MOTION%PINT%HELIUM%RDF%SOLUTE_HE", & 2111 l_val=helium%rdf_sol_he) 2112 CALL section_vals_val_get( & 2113 helium%input, & 2114 "MOTION%PINT%HELIUM%RDF%HE_HE", & 2115 l_val=helium%rdf_he_he) 2116 2117 ! Prevent He_He Rdfs for a single he atom: 2118 IF (helium%atoms <= 1) THEN 2119 helium%rdf_he_he = .FALSE. 2120 END IF 2121 2122 helium%rdf_num = 0 2123 IF (helium%rdf_sol_he) THEN 2124 IF (helium%solute_present) THEN 2125 ! get number of centers from solute to store solute positions 2126 ALLOCATE (helium%rdf_centers(helium%beads, helium%solute_atoms*3)) 2127 helium%rdf_centers(:, :) = 0.0_dp 2128 helium%rdf_num = helium%solute_atoms 2129 ELSE 2130 helium%rdf_sol_he = .FALSE. 2131 END IF 2132 END IF 2133 2134 IF (helium%rdf_he_he) THEN 2135 helium%rdf_num = helium%rdf_num + 1 2136 END IF 2137 2138 ! set the flag for RDF and either proceed or return 2139 IF (helium%rdf_num > 0) THEN 2140 helium%rdf_present = .TRUE. 2141 ELSE 2142 helium%rdf_present = .FALSE. 2143 err_str = "HELIUM RDFs requested, but no valid choice of "// & 2144 "parameters specified. No RDF will be computed." 2145 CPWARN(err_str) 2146 RETURN 2147 END IF 2148 2149 ! set the maximum RDF range 2150 CALL section_vals_val_get( & 2151 helium%input, & 2152 "MOTION%PINT%HELIUM%RDF%MAXR", & 2153 explicit=explicit) 2154 IF (explicit) THEN 2155 ! use the value explicitly set in the input 2156 CALL section_vals_val_get( & 2157 helium%input, & 2158 "MOTION%PINT%HELIUM%RDF%MAXR", & 2159 r_val=helium%rdf_maxr) 2160 ELSE 2161 ! use the default value 2162 CALL section_vals_val_get( & 2163 helium%input, & 2164 "MOTION%PINT%HELIUM%DROPLET_RADIUS", & 2165 explicit=explicit) 2166 IF (explicit) THEN 2167 ! use the droplet radius 2168 IF (helium%solute_present .AND. .NOT. helium%periodic) THEN 2169 ! COM of solute is used as center of the box. 2170 ! Therfore distances became larger then droplet_radius 2171 ! when parts of the solute are on opposite site of 2172 ! COM and helium. 2173 ! Use two times droplet_radius for security: 2174 helium%rdf_maxr = helium%droplet_radius*2.0_dp 2175 ELSE 2176 helium%rdf_maxr = helium%droplet_radius 2177 END IF 2178 ELSE 2179 ! use cell_size and cell_shape 2180 ! (they are set regardless of us being periodic or not) 2181 SELECT CASE (helium%cell_shape) 2182 CASE (helium_cell_shape_cube) 2183 helium%rdf_maxr = helium%cell_size/2.0_dp 2184 CASE (helium_cell_shape_octahedron) 2185 helium%rdf_maxr = helium%cell_size*SQRT(3.0_dp)/4.0_dp 2186 CASE DEFAULT 2187 helium%rdf_maxr = 0.0_dp 2188 WRITE (stmp, *) helium%cell_shape 2189 err_str = "cell shape unknown ("//TRIM(ADJUSTL(stmp))//")" 2190 CPABORT(err_str) 2191 END SELECT 2192 END IF 2193 END IF 2194 2195 ! get number of bins and set bin spacing 2196 CALL section_vals_val_get( & 2197 helium%input, & 2198 "MOTION%PINT%HELIUM%RDF%NBIN", & 2199 i_val=helium%rdf_nbin) 2200 helium%rdf_delr = helium%rdf_maxr/REAL(helium%rdf_nbin, dp) 2201 2202 ! get the weighting factor 2203 CALL section_vals_val_get( & 2204 helium%input, & 2205 "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", & 2206 i_val=helium%rdf_iweight) 2207 2208 ! allocate and initialize memory for RDF storage 2209 ii = helium%rdf_num 2210 ij = helium%rdf_nbin 2211 ALLOCATE (helium%rdf_inst(ii, ij)) 2212 ALLOCATE (helium%rdf_accu(ii, ij)) 2213 ALLOCATE (helium%rdf_rstr(ii, ij)) 2214 helium%rdf_inst(:, :) = 0.0_dp 2215 helium%rdf_accu(:, :) = 0.0_dp 2216 helium%rdf_rstr(:, :) = 0.0_dp 2217 2218 RETURN 2219 END SUBROUTINE helium_rdf_init 2220 2221! *************************************************************************** 2222!> \brief Restore the RDFs from the input structure 2223!> \param helium_env ... 2224!> \date 2011-06-22 2225!> \par History 2226!> 2016-07-14 Modified to work with independent helium_env [cschran] 2227!> 2018-10-19 Changed to bead-wise RDFs of solute-He and He-He [cschran] 2228!> \author Lukasz Walewski 2229! ************************************************************************************************** 2230 SUBROUTINE helium_rdf_restore(helium_env) 2231 2232 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 2233 2234 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rdf_restore', & 2235 routineP = moduleN//':'//routineN 2236 2237 CHARACTER(len=default_string_length) :: stmp1, stmp2 2238 CHARACTER(len=max_line_length) :: err_str 2239 INTEGER :: ii, ij, itmp, k, msglen 2240 LOGICAL :: explicit, ltmp 2241 LOGICAL, DIMENSION(:, :), POINTER :: m 2242 REAL(KIND=dp), DIMENSION(:), POINTER :: message 2243 REAL(KIND=dp), DIMENSION(:, :), POINTER :: f 2244 2245 CALL section_vals_val_get(helium_env(1)%helium%input, & 2246 "MOTION%PINT%HELIUM%AVERAGES%RDF", & 2247 explicit=explicit) 2248 IF (explicit) THEN 2249 NULLIFY (message) 2250 CALL section_vals_val_get(helium_env(1)%helium%input, & 2251 "MOTION%PINT%HELIUM%AVERAGES%RDF", & 2252 r_vals=message) 2253 msglen = SIZE(message) 2254 itmp = SIZE(helium_env(1)%helium%rdf_rstr) 2255 ltmp = (msglen .EQ. itmp) 2256 IF (.NOT. ltmp) THEN 2257 stmp1 = "" 2258 WRITE (stmp1, *) msglen 2259 stmp2 = "" 2260 WRITE (stmp2, *) itmp 2261 err_str = "Size of the RDF array in the input ("// & 2262 TRIM(ADJUSTL(stmp1))// & 2263 ") .NE. that in the runtime ("// & 2264 TRIM(ADJUSTL(stmp2))//")." 2265 CPABORT(err_str) 2266 END IF 2267 ELSE 2268 RETURN 2269 END IF 2270 2271 ii = helium_env(1)%helium%rdf_num 2272 ij = helium_env(1)%helium%rdf_nbin 2273 NULLIFY (m, f) 2274 ALLOCATE (m(ii, ij)) 2275 ALLOCATE (f(ii, ij)) 2276 m(:, :) = .TRUE. 2277 f(:, :) = 0.0_dp 2278 2279 DO k = 1, SIZE(helium_env) 2280 helium_env(k)%helium%rdf_rstr(:, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f) 2281 CALL section_vals_val_get(helium_env(k)%helium%input, & 2282 "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", & 2283 i_val=helium_env(k)%helium%rdf_iweight) 2284 END DO 2285 2286 DEALLOCATE (f, m) 2287 NULLIFY (message) 2288 2289 RETURN 2290 END SUBROUTINE helium_rdf_restore 2291 2292! *************************************************************************** 2293!> \brief Release/deallocate RDF related data structures 2294!> \param helium ... 2295!> \date 2014-02-25 2296!> \author Lukasz Walewski 2297! ************************************************************************************************** 2298 SUBROUTINE helium_rdf_release(helium) 2299 2300 TYPE(helium_solvent_type), POINTER :: helium 2301 2302 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rdf_release', & 2303 routineP = moduleN//':'//routineN 2304 2305 IF (helium%rdf_present) THEN 2306 2307 DEALLOCATE ( & 2308 helium%rdf_rstr, & 2309 helium%rdf_accu, & 2310 helium%rdf_inst) 2311 2312 NULLIFY ( & 2313 helium%rdf_rstr, & 2314 helium%rdf_accu, & 2315 helium%rdf_inst) 2316 2317 IF (helium%rdf_sol_he) THEN 2318 DEALLOCATE (helium%rdf_centers) 2319 NULLIFY (helium%rdf_centers) 2320 END IF 2321 2322 END IF 2323 2324 RETURN 2325 END SUBROUTINE helium_rdf_release 2326 2327! *************************************************************************** 2328!> \brief Check whether property <prop> is activated in the input structure 2329!> \param helium ... 2330!> \param prop ... 2331!> \return ... 2332!> \date 2014-06-26 2333!> \author Lukasz Walewski 2334!> \note The property is controlled by two items in the input structure, 2335!> the printkey and the control section. Two settings result in 2336!> the property being considered active: 2337!> 1. printkey is on at the given print level 2338!> 2. control section is explicit and on 2339!> If the property is considered active it should be calculated 2340!> and printed through out the run. 2341! ************************************************************************************************** 2342 FUNCTION helium_property_active(helium, prop) RESULT(is_active) 2343 2344 TYPE(helium_solvent_type), POINTER :: helium 2345 CHARACTER(len=*) :: prop 2346 LOGICAL :: is_active 2347 2348 CHARACTER(len=*), PARAMETER :: routineN = 'helium_property_active', & 2349 routineP = moduleN//':'//routineN 2350 2351 CHARACTER(len=default_string_length) :: input_path 2352 INTEGER :: print_level 2353 LOGICAL :: explicit, is_on 2354 TYPE(cp_logger_type), POINTER :: logger 2355 TYPE(section_vals_type), POINTER :: print_key, section 2356 2357 is_active = .FALSE. 2358 NULLIFY (logger) 2359 logger => cp_get_default_logger() 2360 2361 ! if the printkey is on at this runlevel consider prop to be active 2362 NULLIFY (print_key) 2363 input_path = "MOTION%PINT%HELIUM%PRINT%"//TRIM(ADJUSTL(prop)) 2364 print_key => section_vals_get_subs_vals( & 2365 helium%input, & 2366 input_path) 2367 is_on = cp_printkey_is_on( & 2368 iteration_info=logger%iter_info, & 2369 print_key=print_key) 2370 IF (is_on) THEN 2371 is_active = .TRUE. 2372 RETURN 2373 END IF 2374 2375 ! if the control section is explicit and on consider prop to be active 2376 ! and activate the printkey 2377 is_active = .FALSE. 2378 NULLIFY (section) 2379 input_path = "MOTION%PINT%HELIUM%"//TRIM(ADJUSTL(prop)) 2380 section => section_vals_get_subs_vals( & 2381 helium%input, & 2382 input_path) 2383 CALL section_vals_get(section, explicit=explicit) 2384 IF (explicit) THEN 2385 ! control section explicitly present, continue checking 2386 CALL section_vals_val_get( & 2387 section, & 2388 "_SECTION_PARAMETERS_", & 2389 l_val=is_on) 2390 IF (is_on) THEN 2391 ! control section is explicit and on, activate the property 2392 is_active = .TRUE. 2393 ! activate the corresponding print_level as well 2394 print_level = logger%iter_info%print_level 2395 CALL section_vals_val_set( & 2396 print_key, & 2397 "_SECTION_PARAMETERS_", & 2398 i_val=print_level) 2399 END IF 2400 END IF 2401 2402 RETURN 2403 END FUNCTION helium_property_active 2404 2405! *************************************************************************** 2406!> \brief Create the density related data structures and initialize 2407!> \param helium ... 2408!> \date 2014-02-25 2409!> \author Lukasz Walewski 2410! ************************************************************************************************** 2411 SUBROUTINE helium_rho_property_init(helium) 2412 2413 TYPE(helium_solvent_type), POINTER :: helium 2414 2415 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rho_property_init', & 2416 routineP = moduleN//':'//routineN 2417 2418 INTEGER :: nc 2419 2420 ALLOCATE (helium%rho_property(rho_num)) 2421 2422 helium%rho_property(rho_atom_number)%name = 'Atom number density' 2423 nc = 1 2424 helium%rho_property(rho_atom_number)%num_components = nc 2425 ALLOCATE (helium%rho_property(rho_atom_number)%filename_suffix(nc)) 2426 ALLOCATE (helium%rho_property(rho_atom_number)%component_name(nc)) 2427 ALLOCATE (helium%rho_property(rho_atom_number)%component_index(nc)) 2428 helium%rho_property(rho_atom_number)%filename_suffix(1) = 'an' 2429 helium%rho_property(rho_atom_number)%component_name(1) = '' 2430 helium%rho_property(rho_atom_number)%component_index(:) = 0 2431 2432 helium%rho_property(rho_projected_area)%name = 'Projected area squared density, A*A(r)' 2433 nc = 3 2434 helium%rho_property(rho_projected_area)%num_components = nc 2435 ALLOCATE (helium%rho_property(rho_projected_area)%filename_suffix(nc)) 2436 ALLOCATE (helium%rho_property(rho_projected_area)%component_name(nc)) 2437 ALLOCATE (helium%rho_property(rho_projected_area)%component_index(nc)) 2438 helium%rho_property(rho_projected_area)%filename_suffix(1) = 'pa_x' 2439 helium%rho_property(rho_projected_area)%filename_suffix(2) = 'pa_y' 2440 helium%rho_property(rho_projected_area)%filename_suffix(3) = 'pa_z' 2441 helium%rho_property(rho_projected_area)%component_name(1) = 'component x' 2442 helium%rho_property(rho_projected_area)%component_name(2) = 'component y' 2443 helium%rho_property(rho_projected_area)%component_name(3) = 'component z' 2444 helium%rho_property(rho_projected_area)%component_index(:) = 0 2445 2446 helium%rho_property(rho_winding_number)%name = 'Winding number squared density, W*W(r)' 2447 nc = 3 2448 helium%rho_property(rho_winding_number)%num_components = nc 2449 ALLOCATE (helium%rho_property(rho_winding_number)%filename_suffix(nc)) 2450 ALLOCATE (helium%rho_property(rho_winding_number)%component_name(nc)) 2451 ALLOCATE (helium%rho_property(rho_winding_number)%component_index(nc)) 2452 helium%rho_property(rho_winding_number)%filename_suffix(1) = 'wn_x' 2453 helium%rho_property(rho_winding_number)%filename_suffix(2) = 'wn_y' 2454 helium%rho_property(rho_winding_number)%filename_suffix(3) = 'wn_z' 2455 helium%rho_property(rho_winding_number)%component_name(1) = 'component x' 2456 helium%rho_property(rho_winding_number)%component_name(2) = 'component y' 2457 helium%rho_property(rho_winding_number)%component_name(3) = 'component z' 2458 helium%rho_property(rho_winding_number)%component_index(:) = 0 2459 2460 helium%rho_property(rho_winding_cycle)%name = 'Winding number squared density, W^2(r)' 2461 nc = 3 2462 helium%rho_property(rho_winding_cycle)%num_components = nc 2463 ALLOCATE (helium%rho_property(rho_winding_cycle)%filename_suffix(nc)) 2464 ALLOCATE (helium%rho_property(rho_winding_cycle)%component_name(nc)) 2465 ALLOCATE (helium%rho_property(rho_winding_cycle)%component_index(nc)) 2466 helium%rho_property(rho_winding_cycle)%filename_suffix(1) = 'wc_x' 2467 helium%rho_property(rho_winding_cycle)%filename_suffix(2) = 'wc_y' 2468 helium%rho_property(rho_winding_cycle)%filename_suffix(3) = 'wc_z' 2469 helium%rho_property(rho_winding_cycle)%component_name(1) = 'component x' 2470 helium%rho_property(rho_winding_cycle)%component_name(2) = 'component y' 2471 helium%rho_property(rho_winding_cycle)%component_name(3) = 'component z' 2472 helium%rho_property(rho_winding_cycle)%component_index(:) = 0 2473 2474 helium%rho_property(rho_moment_of_inertia)%name = 'Moment of inertia' 2475 nc = 3 2476 helium%rho_property(rho_moment_of_inertia)%num_components = nc 2477 ALLOCATE (helium%rho_property(rho_moment_of_inertia)%filename_suffix(nc)) 2478 ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_name(nc)) 2479 ALLOCATE (helium%rho_property(rho_moment_of_inertia)%component_index(nc)) 2480 helium%rho_property(rho_moment_of_inertia)%filename_suffix(1) = 'mi_x' 2481 helium%rho_property(rho_moment_of_inertia)%filename_suffix(2) = 'mi_y' 2482 helium%rho_property(rho_moment_of_inertia)%filename_suffix(3) = 'mi_z' 2483 helium%rho_property(rho_moment_of_inertia)%component_name(1) = 'component x' 2484 helium%rho_property(rho_moment_of_inertia)%component_name(2) = 'component y' 2485 helium%rho_property(rho_moment_of_inertia)%component_name(3) = 'component z' 2486 helium%rho_property(rho_moment_of_inertia)%component_index(:) = 0 2487 2488 helium%rho_property(:)%is_calculated = .FALSE. 2489 2490 RETURN 2491 END SUBROUTINE helium_rho_property_init 2492 2493! *************************************************************************** 2494!> \brief Create the density related data structures and initialize 2495!> \param helium ... 2496!> \date 2014-02-25 2497!> \author Lukasz Walewski 2498! ************************************************************************************************** 2499 SUBROUTINE helium_rho_init(helium) 2500 2501 TYPE(helium_solvent_type), POINTER :: helium 2502 2503 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rho_init', & 2504 routineP = moduleN//':'//routineN 2505 2506 INTEGER :: ii, itmp, jtmp 2507 LOGICAL :: explicit, ltmp 2508 2509 CALL helium_rho_property_init(helium) 2510 2511 helium%rho_num_act = 0 2512 2513 ! check for atom number density 2514 CALL section_vals_val_get( & 2515 helium%input, & 2516 "MOTION%PINT%HELIUM%RHO%ATOM_NUMBER", & 2517 l_val=ltmp) 2518 IF (ltmp) THEN 2519 helium%rho_property(rho_atom_number)%is_calculated = .TRUE. 2520 helium%rho_num_act = helium%rho_num_act + 1 2521 helium%rho_property(rho_atom_number)%component_index(1) = helium%rho_num_act 2522 END IF 2523 2524 ! check for projected area density 2525 CALL section_vals_val_get( & 2526 helium%input, & 2527 "MOTION%PINT%HELIUM%RHO%PROJECTED_AREA_2", & 2528 l_val=ltmp) 2529 IF (ltmp) THEN 2530 helium%rho_property(rho_projected_area)%is_calculated = .TRUE. 2531 DO ii = 1, helium%rho_property(rho_projected_area)%num_components 2532 helium%rho_num_act = helium%rho_num_act + 1 2533 helium%rho_property(rho_projected_area)%component_index(ii) = helium%rho_num_act 2534 END DO 2535 END IF 2536 2537 ! check for winding number density 2538 CALL section_vals_val_get( & 2539 helium%input, & 2540 "MOTION%PINT%HELIUM%RHO%WINDING_NUMBER_2", & 2541 l_val=ltmp) 2542 IF (ltmp) THEN 2543 helium%rho_property(rho_winding_number)%is_calculated = .TRUE. 2544 DO ii = 1, helium%rho_property(rho_winding_number)%num_components 2545 helium%rho_num_act = helium%rho_num_act + 1 2546 helium%rho_property(rho_winding_number)%component_index(ii) = helium%rho_num_act 2547 END DO 2548 END IF 2549 2550 ! check for winding cycle density 2551 CALL section_vals_val_get( & 2552 helium%input, & 2553 "MOTION%PINT%HELIUM%RHO%WINDING_CYCLE_2", & 2554 l_val=ltmp) 2555 IF (ltmp) THEN 2556 helium%rho_property(rho_winding_cycle)%is_calculated = .TRUE. 2557 DO ii = 1, helium%rho_property(rho_winding_cycle)%num_components 2558 helium%rho_num_act = helium%rho_num_act + 1 2559 helium%rho_property(rho_winding_cycle)%component_index(ii) = helium%rho_num_act 2560 END DO 2561 END IF 2562 2563 ! check for moment of inertia density 2564 CALL section_vals_val_get( & 2565 helium%input, & 2566 "MOTION%PINT%HELIUM%RHO%MOMENT_OF_INERTIA", & 2567 l_val=ltmp) 2568 IF (ltmp) THEN 2569 helium%rho_property(rho_moment_of_inertia)%is_calculated = .TRUE. 2570 DO ii = 1, helium%rho_property(rho_moment_of_inertia)%num_components 2571 helium%rho_num_act = helium%rho_num_act + 1 2572 helium%rho_property(rho_moment_of_inertia)%component_index(ii) = helium%rho_num_act 2573 END DO 2574 END IF 2575 2576 ! set the cube dimensions, etc (common to all estimators) 2577 helium%rho_maxr = helium%cell_size 2578 CALL section_vals_val_get( & 2579 helium%input, & 2580 "MOTION%PINT%HELIUM%RHO%NBIN", & 2581 i_val=helium%rho_nbin) 2582 helium%rho_delr = helium%rho_maxr/REAL(helium%rho_nbin, dp) 2583 2584 ! check for optional estimators based on winding paths 2585 helium%rho_num_min_len_wdg = 0 2586 CALL section_vals_val_get( & 2587 helium%input, & 2588 "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_WDG", & 2589 explicit=explicit) 2590 IF (explicit) THEN 2591 NULLIFY (helium%rho_min_len_wdg_vals) 2592 CALL section_vals_val_get( & 2593 helium%input, & 2594 "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_WDG", & 2595 i_vals=helium%rho_min_len_wdg_vals) 2596 itmp = SIZE(helium%rho_min_len_wdg_vals) 2597 IF (itmp .GT. 0) THEN 2598 helium%rho_num_min_len_wdg = itmp 2599 helium%rho_num_act = helium%rho_num_act + itmp 2600 END IF 2601 END IF 2602 2603 ! check for optional estimators based on non-winding paths 2604 helium%rho_num_min_len_non = 0 2605 CALL section_vals_val_get( & 2606 helium%input, & 2607 "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_NON", & 2608 explicit=explicit) 2609 IF (explicit) THEN 2610 NULLIFY (helium%rho_min_len_non_vals) 2611 CALL section_vals_val_get( & 2612 helium%input, & 2613 "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_NON", & 2614 i_vals=helium%rho_min_len_non_vals) 2615 itmp = SIZE(helium%rho_min_len_non_vals) 2616 IF (itmp .GT. 0) THEN 2617 helium%rho_num_min_len_non = itmp 2618 helium%rho_num_act = helium%rho_num_act + itmp 2619 END IF 2620 END IF 2621 2622 ! check for optional estimators based on all paths 2623 helium%rho_num_min_len_all = 0 2624 CALL section_vals_val_get( & 2625 helium%input, & 2626 "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_ALL", & 2627 explicit=explicit) 2628 IF (explicit) THEN 2629 NULLIFY (helium%rho_min_len_all_vals) 2630 CALL section_vals_val_get( & 2631 helium%input, & 2632 "MOTION%PINT%HELIUM%RHO%MIN_CYCLE_LENGTHS_ALL", & 2633 i_vals=helium%rho_min_len_all_vals) 2634 itmp = SIZE(helium%rho_min_len_all_vals) 2635 IF (itmp .GT. 0) THEN 2636 helium%rho_num_min_len_all = itmp 2637 helium%rho_num_act = helium%rho_num_act + itmp 2638 END IF 2639 END IF 2640 2641 ! get the weighting factor 2642 CALL section_vals_val_get( & 2643 helium%input, & 2644 "MOTION%PINT%HELIUM%AVERAGES%IWEIGHT", & 2645 i_val=helium%rho_iweight) 2646 2647 ! allocate and initialize memory for density storage 2648 itmp = helium%rho_nbin 2649 jtmp = helium%rho_num_act 2650 ALLOCATE (helium%rho_inst(jtmp, itmp, itmp, itmp)) 2651 ALLOCATE (helium%rho_accu(jtmp, itmp, itmp, itmp)) 2652 ALLOCATE (helium%rho_rstr(jtmp, itmp, itmp, itmp)) 2653 ALLOCATE (helium%rho_incr(jtmp, helium%atoms, helium%beads)) 2654 2655 helium%rho_incr(:, :, :) = 0.0_dp 2656 helium%rho_inst(:, :, :, :) = 0.0_dp 2657 helium%rho_accu(:, :, :, :) = 0.0_dp 2658 helium%rho_rstr(:, :, :, :) = 0.0_dp 2659 2660 RETURN 2661 END SUBROUTINE helium_rho_init 2662 2663! *************************************************************************** 2664!> \brief Restore the densities from the input structure. 2665!> \param helium_env ... 2666!> \date 2011-06-22 2667!> \par History 2668!> 2016-07-14 Modified to work with independent helium_env [cschran] 2669!> \author Lukasz Walewski 2670! ************************************************************************************************** 2671 SUBROUTINE helium_rho_restore(helium_env) 2672 2673 TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env 2674 2675 CHARACTER(len=*), PARAMETER :: routineN = 'helium_rho_restore', & 2676 routineP = moduleN//':'//routineN 2677 2678 CHARACTER(len=default_string_length) :: stmp1, stmp2 2679 CHARACTER(len=max_line_length) :: err_str 2680 INTEGER :: itmp, k, msglen 2681 LOGICAL :: explicit, ltmp 2682 LOGICAL, DIMENSION(:, :, :, :), POINTER :: m 2683 REAL(KIND=dp), DIMENSION(:), POINTER :: message 2684 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: f 2685 2686 CALL section_vals_val_get(helium_env(1)%helium%input, & 2687 "MOTION%PINT%HELIUM%AVERAGES%RHO", & 2688 explicit=explicit) 2689 IF (explicit) THEN 2690 NULLIFY (message) 2691 CALL section_vals_val_get(helium_env(1)%helium%input, & 2692 "MOTION%PINT%HELIUM%AVERAGES%RHO", & 2693 r_vals=message) 2694 msglen = SIZE(message) 2695 itmp = SIZE(helium_env(1)%helium%rho_rstr) 2696 ltmp = (msglen .EQ. itmp) 2697 IF (.NOT. ltmp) THEN 2698 stmp1 = "" 2699 WRITE (stmp1, *) msglen 2700 stmp2 = "" 2701 WRITE (stmp2, *) itmp 2702 err_str = "Size of the S density array in the input ("// & 2703 TRIM(ADJUSTL(stmp1))// & 2704 ") .NE. that in the runtime ("// & 2705 TRIM(ADJUSTL(stmp2))//")." 2706 CPABORT(err_str) 2707 END IF 2708 ELSE 2709 RETURN 2710 END IF 2711 2712 itmp = helium_env(1)%helium%rho_nbin 2713 NULLIFY (m, f) 2714 ALLOCATE (m(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp)) 2715 ALLOCATE (f(helium_env(1)%helium%rho_num_act, itmp, itmp, itmp)) 2716 m(:, :, :, :) = .TRUE. 2717 f(:, :, :, :) = 0.0_dp 2718 2719 DO k = 1, SIZE(helium_env) 2720 helium_env(k)%helium%rho_rstr(:, :, :, :) = UNPACK(message(1:msglen), MASK=m, FIELD=f) 2721 END DO 2722 2723 DEALLOCATE (f, m) 2724 NULLIFY (message) 2725 2726 RETURN 2727 END SUBROUTINE helium_rho_restore 2728 2729! *************************************************************************** 2730!> \brief Count atoms of different types and store their global indices. 2731!> \param helium ... 2732!> \param pint_env ... 2733!> \author Lukasz Walewski 2734!> \note Arrays ALLOCATEd here are (should be) DEALLOCATEd in 2735!> helium_release. 2736! ************************************************************************************************** 2737 SUBROUTINE helium_set_solute_indices(helium, pint_env) 2738 TYPE(helium_solvent_type), POINTER :: helium 2739 TYPE(pint_env_type), POINTER :: pint_env 2740 2741 CHARACTER(len=*), PARAMETER :: routineN = 'helium_set_solute_indices', & 2742 routineP = moduleN//':'//routineN 2743 2744 INTEGER :: iatom, natoms 2745 REAL(KIND=dp) :: mass 2746 TYPE(cp_subsys_type), POINTER :: my_subsys 2747 TYPE(f_env_type), POINTER :: my_f_env 2748 TYPE(particle_list_type), POINTER :: my_particles 2749 2750! set up my_particles structure 2751 2752 NULLIFY (my_f_env, my_subsys, my_particles) 2753 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, & 2754 f_env=my_f_env) 2755 CALL force_env_get(force_env=my_f_env%force_env, subsys=my_subsys) 2756 CALL cp_subsys_get(my_subsys, particles=my_particles) 2757 CALL f_env_rm_defaults(my_f_env) 2758 2759 natoms = helium%solute_atoms 2760 NULLIFY (helium%solute_element) 2761 ALLOCATE (helium%solute_element(natoms)) 2762 2763 ! find out how many different atomic types are there 2764 helium%enum = 0 2765 DO iatom = 1, natoms 2766 CALL get_atomic_kind(my_particles%els(iatom)%atomic_kind, & 2767 mass=mass, & 2768 element_symbol=helium%solute_element(iatom)) 2769 END DO 2770 2771 RETURN 2772 END SUBROUTINE helium_set_solute_indices 2773 2774! *************************************************************************** 2775!> \brief Sets helium%solute_cell based on the solute's force_env. 2776!> \param helium ... 2777!> \param pint_env ... 2778!> \author Lukasz Walewski 2779!> \note The simulation cell for the solvated molecule is taken from force_env 2780!> which should assure that we get proper cell dimensions regardless of 2781!> the method used for the solute (QS, FIST). Helium solvent needs the 2782!> solute's cell dimensions to calculate the solute-solvent distances 2783!> correctly. 2784!> \note At the moment only orthorhombic cells are supported. 2785! ************************************************************************************************** 2786 SUBROUTINE helium_set_solute_cell(helium, pint_env) 2787 TYPE(helium_solvent_type), POINTER :: helium 2788 TYPE(pint_env_type), POINTER :: pint_env 2789 2790 CHARACTER(len=*), PARAMETER :: routineN = 'helium_set_solute_cell', & 2791 routineP = moduleN//':'//routineN 2792 2793 LOGICAL :: my_orthorhombic 2794 TYPE(cell_type), POINTER :: my_cell 2795 TYPE(f_env_type), POINTER :: my_f_env 2796 2797! get the cell structure from pint_env 2798 2799 NULLIFY (my_f_env, my_cell) 2800 CALL f_env_add_defaults(f_env_id=pint_env%replicas%f_env_id, & 2801 f_env=my_f_env) 2802 CALL force_env_get(force_env=my_f_env%force_env, cell=my_cell) 2803 CALL f_env_rm_defaults(my_f_env) 2804 2805 CALL get_cell(my_cell, orthorhombic=my_orthorhombic) 2806 IF (.NOT. my_orthorhombic) THEN 2807 CPABORT("Helium solvent not implemented for non-orthorhombic cells.") 2808 ELSE 2809 helium%solute_cell => my_cell 2810 END IF 2811 2812 RETURN 2813 END SUBROUTINE helium_set_solute_cell 2814 2815END MODULE helium_methods 2816