1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6MODULE input_restart_force_eval 7 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 8 USE atomic_kind_types, ONLY: get_atomic_kind 9 USE cell_types, ONLY: cell_type,& 10 real_to_scaled 11 USE cp_control_types, ONLY: dft_control_type 12 USE cp_files, ONLY: close_file,& 13 open_file 14 USE cp_linked_list_input, ONLY: cp_sll_val_create,& 15 cp_sll_val_get_length,& 16 cp_sll_val_type 17 USE cp_para_types, ONLY: cp_para_env_type 18 USE cp_subsys_types, ONLY: cp_subsys_get,& 19 cp_subsys_type 20 USE cp_units, ONLY: cp_unit_from_cp2k 21 USE distribution_1d_types, ONLY: distribution_1d_type 22 USE eip_environment_types, ONLY: eip_env_get 23 USE force_env_types, ONLY: force_env_get,& 24 force_env_type,& 25 multiple_fe_list,& 26 use_eip_force,& 27 use_qmmm,& 28 use_qs_force 29 USE input_constants, ONLY: do_coord_cp2k,& 30 ehrenfest,& 31 mol_dyn_run,& 32 mon_car_run,& 33 pint_run,& 34 use_rt_restart 35 USE input_cp2k_restarts_util, ONLY: section_velocity_val_set 36 USE input_restart_rng, ONLY: section_rng_val_set 37 USE input_section_types, ONLY: & 38 section_get_keyword_index, section_type, section_vals_add_values, section_vals_get, & 39 section_vals_get_subs_vals, section_vals_get_subs_vals3, section_vals_remove_values, & 40 section_vals_type, section_vals_val_get, section_vals_val_set, section_vals_val_unset 41 USE input_val_types, ONLY: val_create,& 42 val_release,& 43 val_type 44 USE kinds, ONLY: default_string_length,& 45 dp 46 USE message_passing, ONLY: mp_sum 47 USE molecule_kind_types, ONLY: get_molecule_kind 48 USE molecule_list_types, ONLY: molecule_list_type 49 USE molecule_types, ONLY: get_molecule,& 50 molecule_type 51 USE multipole_types, ONLY: multipole_type 52 USE parallel_rng_types, ONLY: dump_rng_stream,& 53 rng_record_length 54 USE particle_list_types, ONLY: particle_list_type 55 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 56 USE qs_environment_types, ONLY: get_qs_env 57 USE string_utilities, ONLY: string_to_ascii 58 USE virial_types, ONLY: virial_type 59#include "./base/base_uses.f90" 60 61 IMPLICIT NONE 62 63 PRIVATE 64 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_restart_force_eval' 66 67 PUBLIC :: update_force_eval, update_subsys 68 69CONTAINS 70 71! ************************************************************************************************** 72!> \brief Updates the force_eval section of the input file 73!> \param force_env ... 74!> \param root_section ... 75!> \param write_binary_restart_file ... 76!> \param respa ... 77!> \par History 78!> 01.2006 created [teo] 79!> \author Teodoro Laino 80! ************************************************************************************************** 81 SUBROUTINE update_force_eval(force_env, root_section, & 82 write_binary_restart_file, respa) 83 TYPE(force_env_type), POINTER :: force_env 84 TYPE(section_vals_type), POINTER :: root_section 85 LOGICAL, INTENT(IN) :: write_binary_restart_file 86 LOGICAL, OPTIONAL :: respa 87 88 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_force_eval', & 89 routineP = moduleN//':'//routineN 90 91 INTEGER :: i_subsys, iforce_eval, myid, nforce_eval 92 INTEGER, DIMENSION(:), POINTER :: i_force_eval 93 LOGICAL :: multiple_subsys, skip_vel_section 94 REAL(KIND=dp), DIMENSION(:), POINTER :: work 95 TYPE(cell_type), POINTER :: cell 96 TYPE(cp_subsys_type), POINTER :: subsys 97 TYPE(dft_control_type), POINTER :: dft_control 98 TYPE(section_vals_type), POINTER :: cell_section, dft_section, & 99 force_env_sections, qmmm_section, & 100 rng_section, subsys_section 101 TYPE(virial_type), POINTER :: virial 102 103 NULLIFY (rng_section, subsys_section, cell_section, virial, subsys, cell, dft_control, work) 104 ! If it's not a dynamical run don't update the velocity section 105 CALL section_vals_val_get(root_section, "GLOBAL%RUN_TYPE", i_val=myid) 106 skip_vel_section = ((myid /= mol_dyn_run) .AND. & 107 (myid /= mon_car_run) .AND. & 108 (myid /= pint_run) .AND. & 109 (myid /= ehrenfest)) 110 111 ! Go on updatig the force_env_section 112 force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL") 113 CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, nforce_eval) 114 ! The update of the input MUST be realized only on the main force_eval 115 ! All the others will be left not updated because there is no real need to update them... 116 iforce_eval = 1 117 subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", & 118 i_rep_section=i_force_eval(iforce_eval)) 119 CALL update_subsys(subsys_section, force_env, skip_vel_section, & 120 write_binary_restart_file) 121 122 ! For RESPA we need to update all subsystems 123 IF (PRESENT(respa)) THEN 124 IF (respa) THEN 125 DO i_subsys = 1, 2 126 iforce_eval = i_subsys 127 subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", & 128 i_rep_section=i_force_eval(iforce_eval)) 129 CALL update_subsys(subsys_section, force_env, skip_vel_section, & 130 write_binary_restart_file) 131 ENDDO 132 ENDIF 133 END IF 134 135 rng_section => section_vals_get_subs_vals(subsys_section, "RNG_INIT") 136 CALL update_rng_particle(rng_section, force_env) 137 138 qmmm_section => section_vals_get_subs_vals3(force_env_sections, "QMMM", & 139 i_rep_section=i_force_eval(iforce_eval)) 140 CALL update_qmmm(qmmm_section, force_env) 141 142 ! Non-mixed CDFT calculation: update constraint Lagrangian and counter 143 IF (nforce_eval == 1) THEN 144 IF (ASSOCIATED(force_env%qs_env)) THEN 145 CALL get_qs_env(force_env%qs_env, dft_control=dft_control) 146 IF (dft_control%qs_control%cdft) THEN 147 dft_section => section_vals_get_subs_vals3(force_env_sections, "DFT", & 148 i_rep_section=i_force_eval(iforce_eval)) 149 CALL section_vals_val_set(dft_section, "QS%CDFT%COUNTER", & 150 i_val=dft_control%qs_control%cdft_control%ienergy) 151 ALLOCATE (work(SIZE(dft_control%qs_control%cdft_control%strength))) 152 work = dft_control%qs_control%cdft_control%strength 153 CALL section_vals_val_set(dft_section, "QS%CDFT%STRENGTH", & 154 r_vals_ptr=work) 155 END IF 156 END IF 157 END IF 158 ! And now update only the cells of all other force_eval 159 ! This is to make consistent for cell variable runs.. 160 IF (nforce_eval > 1) THEN 161 CALL force_env_get(force_env, subsys=subsys, cell=cell) 162 CALL cp_subsys_get(subsys, virial=virial) 163 CALL section_vals_val_get(root_section, "MULTIPLE_FORCE_EVALS%MULTIPLE_SUBSYS", l_val=multiple_subsys) 164 IF (virial%pv_availability .AND. multiple_subsys) THEN 165 DO iforce_eval = 2, nforce_eval 166 subsys_section => section_vals_get_subs_vals3(force_env_sections, "SUBSYS", & 167 i_rep_section=i_force_eval(iforce_eval)) 168 cell_section => section_vals_get_subs_vals(subsys_section, "CELL") 169 CALL update_cell_section(cell, cell_section) 170 END DO 171 END IF 172 ! With mixed CDFT, update value of constraint Lagrangian 173 IF (ASSOCIATED(force_env%mixed_env)) THEN 174 IF (force_env%mixed_env%do_mixed_cdft) THEN 175 DO iforce_eval = 2, nforce_eval 176 dft_section => section_vals_get_subs_vals3(force_env_sections, "DFT", & 177 i_rep_section=i_force_eval(iforce_eval)) 178 ALLOCATE (work(SIZE(force_env%mixed_env%strength, 2))) 179 work = force_env%mixed_env%strength(iforce_eval - 1, :) 180 CALL section_vals_val_set(dft_section, "QS%CDFT%STRENGTH", & 181 r_vals_ptr=work) 182 CALL section_vals_val_set(dft_section, "QS%CDFT%COUNTER", & 183 i_val=force_env%mixed_env%cdft_control%sim_step) 184 END DO 185 END IF 186 END IF 187 END IF 188 189 IF (myid == ehrenfest) CALL section_vals_val_set(root_section, "FORCE_EVAL%DFT%REAL_TIME_PROPAGATION%INITIAL_WFN", & 190 i_val=use_rt_restart) 191 DEALLOCATE (i_force_eval) 192 193 END SUBROUTINE update_force_eval 194 195! ************************************************************************************************** 196!> \brief Updates the qmmm section if needed 197!> \param qmmm_section ... 198!> \param force_env ... 199!> \par History 200!> 08.2007 created [teo] 201!> \author Teodoro Laino 202! ************************************************************************************************** 203 SUBROUTINE update_qmmm(qmmm_section, force_env) 204 TYPE(section_vals_type), POINTER :: qmmm_section 205 TYPE(force_env_type), POINTER :: force_env 206 207 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_qmmm', routineP = moduleN//':'//routineN 208 209 LOGICAL :: explicit 210 REAL(KIND=dp), DIMENSION(:), POINTER :: work 211 212 SELECT CASE (force_env%in_use) 213 CASE (use_qmmm) 214 CALL section_vals_get(qmmm_section, explicit=explicit) 215 CPASSERT(explicit) 216 217 ALLOCATE (work(3)) 218 work = force_env%qmmm_env%qm%transl_v 219 CALL section_vals_val_set(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals_ptr=work) 220 END SELECT 221 222 END SUBROUTINE update_qmmm 223 224! ************************************************************************************************** 225!> \brief Updates the rng section of the input file 226!> Write current status of the parallel random number generator (RNG) 227!> \param rng_section ... 228!> \param force_env ... 229!> \par History 230!> 01.2006 created [teo] 231!> \author Teodoro Laino 232! ************************************************************************************************** 233 SUBROUTINE update_rng_particle(rng_section, force_env) 234 235 TYPE(section_vals_type), POINTER :: rng_section 236 TYPE(force_env_type), POINTER :: force_env 237 238 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_rng_particle', & 239 routineP = moduleN//':'//routineN 240 241 CHARACTER(LEN=rng_record_length) :: rng_record 242 INTEGER :: iparticle, iparticle_kind, & 243 iparticle_local, nparticle, & 244 nparticle_kind, nparticle_local 245 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ascii 246 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 247 TYPE(cp_para_env_type), POINTER :: para_env 248 TYPE(cp_subsys_type), POINTER :: subsys 249 TYPE(distribution_1d_type), POINTER :: local_particles 250 TYPE(particle_list_type), POINTER :: particles 251 252 CALL force_env_get(force_env, subsys=subsys, para_env=para_env) 253 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 254 particles=particles) 255 256 IF (ASSOCIATED(local_particles%local_particle_set)) THEN 257 nparticle_kind = atomic_kinds%n_els 258 nparticle = particles%n_els 259 260 ALLOCATE (ascii(rng_record_length, nparticle)) 261 ascii = 0 262 263 DO iparticle = 1, nparticle 264 DO iparticle_kind = 1, nparticle_kind 265 nparticle_local = local_particles%n_el(iparticle_kind) 266 DO iparticle_local = 1, nparticle_local 267 IF (iparticle == local_particles%list(iparticle_kind)%array(iparticle_local)) THEN 268 CALL dump_rng_stream(rng_stream=local_particles%local_particle_set(iparticle_kind)% & 269 rng(iparticle_local)%stream, & 270 rng_record=rng_record) 271 CALL string_to_ascii(rng_record, ascii(:, iparticle)) 272 END IF 273 END DO 274 END DO 275 END DO 276 277 CALL mp_sum(ascii, para_env%group) 278 279 CALL section_rng_val_set(rng_section=rng_section, nsize=nparticle, ascii=ascii) 280 281 DEALLOCATE (ascii) 282 END IF 283 284 END SUBROUTINE update_rng_particle 285 286! ************************************************************************************************** 287!> \brief Updates the subsys section of the input file 288!> \param subsys_section ... 289!> \param force_env ... 290!> \param skip_vel_section ... 291!> \param write_binary_restart_file ... 292!> \par History 293!> 01.2006 created [teo] 294!> \author Teodoro Laino 295! ************************************************************************************************** 296 SUBROUTINE update_subsys(subsys_section, force_env, skip_vel_section, & 297 write_binary_restart_file) 298 TYPE(section_vals_type), POINTER :: subsys_section 299 TYPE(force_env_type), POINTER :: force_env 300 LOGICAL, INTENT(IN) :: skip_vel_section, & 301 write_binary_restart_file 302 303 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_subsys', routineP = moduleN//':'//routineN 304 305 CHARACTER(LEN=default_string_length) :: coord_file_name, unit_str 306 INTEGER :: coord_file_format, handle, output_unit 307 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell 308 LOGICAL :: scale, use_ref_cell 309 REAL(KIND=dp) :: conv_factor 310 TYPE(cell_type), POINTER :: cell 311 TYPE(cp_para_env_type), POINTER :: para_env 312 TYPE(cp_subsys_type), POINTER :: subsys 313 TYPE(molecule_list_type), POINTER :: molecules 314 TYPE(multipole_type), POINTER :: multipoles 315 TYPE(particle_list_type), POINTER :: core_particles, particles, & 316 shell_particles 317 TYPE(section_vals_type), POINTER :: work_section 318 319 NULLIFY (work_section, core_particles, particles, shell_particles, & 320 subsys, cell, para_env, multipoles) 321 CALL timeset(routineN, handle) 322 CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env) 323 324 CALL cp_subsys_get(subsys, particles=particles, molecules=molecules, & 325 shell_particles=shell_particles, core_particles=core_particles, & 326 multipoles=multipoles) 327 328 ! Remove the multiple_unit_cell from the input structure.. at this point we have 329 ! already all the information we need.. 330 ALLOCATE (multiple_unit_cell(3)) 331 multiple_unit_cell = 1 332 CALL section_vals_val_set(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", & 333 i_vals_ptr=multiple_unit_cell) 334 335 ! Coordinates and Velocities 336 work_section => section_vals_get_subs_vals(subsys_section, "COORD") 337 CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str) 338 CALL section_vals_val_get(work_section, "SCALED", l_val=scale) 339 conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 340 IF (.NOT. write_binary_restart_file) THEN 341 CALL section_vals_val_get(subsys_section, "TOPOLOGY%COORD_FILE_FORMAT", & 342 i_val=coord_file_format) 343 IF (coord_file_format == do_coord_cp2k) THEN 344 CALL section_vals_val_get(subsys_section, "TOPOLOGY%COORD_FILE_NAME", & 345 c_val=coord_file_name) 346 output_unit = 0 347 IF (para_env%ionode) THEN 348 CALL open_file(file_name=TRIM(ADJUSTL(coord_file_name)), & 349 file_action="READWRITE", & 350 file_form="FORMATTED", & 351 file_position="REWIND", & 352 file_status="UNKNOWN", & 353 unit_number=output_unit) 354 CALL dump_coordinates_cp2k(particles, molecules, cell, conv_factor, & 355 output_unit=output_unit, & 356 core_or_shell=.FALSE., & 357 scaled_coordinates=scale) 358 CALL close_file(unit_number=output_unit) 359 END IF 360 ELSE 361 CALL section_coord_val_set(work_section, particles, molecules, conv_factor, scale, & 362 cell) 363 END IF 364 END IF 365 CALL section_vals_val_set(subsys_section, "TOPOLOGY%NUMBER_OF_ATOMS", i_val=particles%n_els) 366 work_section => section_vals_get_subs_vals(subsys_section, "VELOCITY") 367 IF (.NOT. skip_vel_section) THEN 368 IF (.NOT. write_binary_restart_file) THEN 369 CALL section_velocity_val_set(work_section, particles, conv_factor=1.0_dp) 370 END IF 371 ELSE 372 CALL section_vals_remove_values(work_section) 373 END IF 374 375 ! Write restart input for core-shell model 376 IF (.NOT. write_binary_restart_file) THEN 377 IF (ASSOCIATED(shell_particles)) THEN 378 work_section => section_vals_get_subs_vals(subsys_section, "SHELL_COORD") 379 CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str) 380 CALL section_vals_val_get(work_section, "SCALED", l_val=scale) 381 conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 382 CALL section_coord_val_set(work_section, shell_particles, molecules, & 383 conv_factor, scale, cell, shell=.TRUE.) 384 IF (.NOT. skip_vel_section) THEN 385 work_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY") 386 CALL section_velocity_val_set(work_section, shell_particles, conv_factor=1.0_dp) 387 END IF 388 END IF 389 IF (ASSOCIATED(core_particles)) THEN 390 work_section => section_vals_get_subs_vals(subsys_section, "CORE_COORD") 391 CALL section_vals_val_get(work_section, "UNIT", c_val=unit_str) 392 CALL section_vals_val_get(work_section, "SCALED", l_val=scale) 393 conv_factor = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 394 CALL section_coord_val_set(work_section, core_particles, molecules, & 395 conv_factor, scale, cell, shell=.TRUE.) 396 IF (.NOT. skip_vel_section) THEN 397 work_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY") 398 CALL section_velocity_val_set(work_section, core_particles, conv_factor=1.0_dp) 399 END IF 400 END IF 401 END IF 402 403 ! Updating cell info 404 CALL force_env_get(force_env, cell=cell) 405 work_section => section_vals_get_subs_vals(subsys_section, "CELL") 406 CALL update_cell_section(cell, cell_section=work_section) 407 408 ! Updating cell_ref info 409 use_ref_cell = .FALSE. 410 SELECT CASE (force_env%in_use) 411 CASE (use_qs_force) 412 CALL get_qs_env(force_env%qs_env, cell_ref=cell, use_ref_cell=use_ref_cell) 413 CASE (use_eip_force) 414 CALL eip_env_get(force_env%eip_env, cell_ref=cell, use_ref_cell=use_ref_cell) 415 END SELECT 416 IF (use_ref_cell) THEN 417 work_section => section_vals_get_subs_vals(subsys_section, "CELL%CELL_REF") 418 CALL update_cell_section(cell, cell_section=work_section) 419 ENDIF 420 421 ! Updating multipoles 422 IF (ASSOCIATED(multipoles)) THEN 423 work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES") 424 DO 425 IF (SIZE(work_section%values, 2) == 1) EXIT 426 CALL section_vals_add_values(work_section) 427 END DO 428 IF (ASSOCIATED(multipoles%dipoles)) THEN 429 work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%DIPOLES") 430 CALL update_dipoles_section(multipoles%dipoles, work_section) 431 END IF 432 IF (ASSOCIATED(multipoles%quadrupoles)) THEN 433 work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%QUADRUPOLES") 434 CALL update_quadrupoles_section(multipoles%quadrupoles, work_section) 435 END IF 436 END IF 437 CALL timestop(handle) 438 END SUBROUTINE update_subsys 439 440! ************************************************************************************************** 441!> \brief Routine to update a cell section 442!> \param cell ... 443!> \param cell_section ... 444!> \author Ole Schuett 445! ************************************************************************************************** 446 SUBROUTINE update_cell_section(cell, cell_section) 447 TYPE(cell_type), POINTER :: cell 448 TYPE(section_vals_type), POINTER :: cell_section 449 450 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_cell_section', & 451 routineP = moduleN//':'//routineN 452 453 INTEGER :: handle 454 INTEGER, DIMENSION(:), POINTER :: iwork 455 REAL(KIND=dp), DIMENSION(:), POINTER :: work 456 457 NULLIFY (work, iwork) 458 CALL timeset(routineN, handle) 459 460 ! CELL VECTORS - A 461 ALLOCATE (work(3)) 462 work(1:3) = cell%hmat(1:3, 1) 463 CALL section_vals_val_set(cell_section, "A", r_vals_ptr=work) 464 465 ! CELL VECTORS - B 466 ALLOCATE (work(3)) 467 work(1:3) = cell%hmat(1:3, 2) 468 CALL section_vals_val_set(cell_section, "B", r_vals_ptr=work) 469 470 ! CELL VECTORS - C 471 ALLOCATE (work(3)) 472 work(1:3) = cell%hmat(1:3, 3) 473 CALL section_vals_val_set(cell_section, "C", r_vals_ptr=work) 474 475 ! MULTIPLE_UNIT_CELL 476 ALLOCATE (iwork(3)) 477 iwork = 1 478 CALL section_vals_val_set(cell_section, "MULTIPLE_UNIT_CELL", i_vals_ptr=iwork) 479 480 ! Unset unused or misleading information 481 CALL section_vals_val_unset(cell_section, "ABC") 482 CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA") 483 484 CALL timestop(handle) 485 END SUBROUTINE update_cell_section 486 487! ************************************************************************************************** 488!> \brief routine to dump coordinates.. fast implementation 489!> \param coord_section ... 490!> \param particles ... 491!> \param molecules ... 492!> \param conv_factor ... 493!> \param scale ... 494!> \param cell ... 495!> \param shell ... 496!> \par History 497!> 02.2006 created [teo] 498!> \author Teodoro Laino 499! ************************************************************************************************** 500 SUBROUTINE section_coord_val_set(coord_section, particles, molecules, conv_factor, & 501 scale, cell, shell) 502 TYPE(section_vals_type), POINTER :: coord_section 503 TYPE(particle_list_type), POINTER :: particles 504 TYPE(molecule_list_type), POINTER :: molecules 505 REAL(KIND=dp) :: conv_factor 506 LOGICAL, INTENT(IN) :: scale 507 TYPE(cell_type), POINTER :: cell 508 LOGICAL, INTENT(IN), OPTIONAL :: shell 509 510 CHARACTER(LEN=*), PARAMETER :: routineN = 'section_coord_val_set', & 511 routineP = moduleN//':'//routineN 512 513 CHARACTER(LEN=2*default_string_length) :: line 514 CHARACTER(LEN=default_string_length) :: my_tag, name 515 INTEGER :: handle, ik, imol, irk, last_atom, Nlist 516 LOGICAL :: ldum, molname_generated, my_shell 517 REAL(KIND=dp), DIMENSION(3) :: r, s 518 TYPE(cp_sll_val_type), POINTER :: new_pos, vals 519 TYPE(molecule_type), POINTER :: molecule_now 520 TYPE(section_type), POINTER :: section 521 TYPE(val_type), POINTER :: my_val, old_val 522 523 CALL timeset(routineN, handle) 524 525 NULLIFY (my_val, old_val, section, vals) 526 my_shell = .FALSE. 527 IF (PRESENT(shell)) my_shell = shell 528 CPASSERT(ASSOCIATED(coord_section)) 529 CPASSERT(coord_section%ref_count > 0) 530 section => coord_section%section 531 ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_") 532 IF (ik == -2) & 533 CALL cp_abort(__LOCATION__, & 534 "section "//TRIM(section%name)//" does not contain keyword "// & 535 "_DEFAULT_KEYWORD_") 536 537 DO 538 IF (SIZE(coord_section%values, 2) == 1) EXIT 539 CALL section_vals_add_values(coord_section) 540 END DO 541 vals => coord_section%values(ik, 1)%list 542 Nlist = 0 543 IF (ASSOCIATED(vals)) THEN 544 Nlist = cp_sll_val_get_length(vals) 545 END IF 546 imol = 0 547 last_atom = 0 548 DO irk = 1, particles%n_els 549 CALL get_atomic_kind(particles%els(irk)%atomic_kind, name=name) 550 IF (my_shell) THEN 551 s = particles%els(irk)%r 552 IF (scale) THEN 553 r = s 554 CALL real_to_scaled(s, r, cell) 555 ELSE 556 s = s*conv_factor 557 END IF 558 WRITE (UNIT=line, FMT="(A,3(1X,ES25.16),1X,I0)") & 559 TRIM(name), s(1:3), particles%els(irk)%atom_index 560 CALL val_create(my_val, lc_val=line) 561 IF (Nlist /= 0) THEN 562 IF (irk == 1) THEN 563 new_pos => vals 564 ELSE 565 new_pos => new_pos%rest 566 END IF 567 old_val => new_pos%first_el 568 CALL val_release(old_val) 569 new_pos%first_el => my_val 570 ELSE 571 IF (irk == 1) THEN 572 NULLIFY (new_pos) 573 CALL cp_sll_val_create(new_pos, first_el=my_val) 574 vals => new_pos 575 ELSE 576 NULLIFY (new_pos%rest) 577 CALL cp_sll_val_create(new_pos%rest, first_el=my_val) 578 new_pos => new_pos%rest 579 END IF 580 END IF 581 NULLIFY (my_val) 582 ELSE 583 IF (last_atom < irk) THEN 584 imol = imol + 1 585 molecule_now => molecules%els(imol) 586 CALL get_molecule(molecule_now, last_atom=last_atom) 587 CALL get_molecule_kind(molecule_now%molecule_kind, molname_generated=molname_generated, & 588 name=my_tag) 589 IF (molname_generated) my_tag = "" 590 END IF 591 ldum = qmmm_ff_precond_only_qm(my_tag) 592 ldum = qmmm_ff_precond_only_qm(name) 593 s = particles%els(irk)%r 594 IF (scale) THEN 595 r = s 596 CALL real_to_scaled(s, r, cell) 597 ELSE 598 s = s*conv_factor 599 END IF 600 IF (LEN_TRIM(my_tag) > 0) THEN 601 WRITE (UNIT=line, FMT="(A,3(1X,ES25.16),1X,A,1X,I0)") & 602 TRIM(name), s(1:3), TRIM(my_tag), imol 603 ELSE 604 WRITE (UNIT=line, FMT="(A,3(1X,ES25.16))") & 605 TRIM(name), s(1:3) 606 END IF 607 CALL val_create(my_val, lc_val=line) 608 609 IF (Nlist /= 0) THEN 610 IF (irk == 1) THEN 611 new_pos => vals 612 ELSE 613 new_pos => new_pos%rest 614 END IF 615 old_val => new_pos%first_el 616 CALL val_release(old_val) 617 new_pos%first_el => my_val 618 ELSE 619 IF (irk == 1) THEN 620 NULLIFY (new_pos) 621 CALL cp_sll_val_create(new_pos, first_el=my_val) 622 vals => new_pos 623 ELSE 624 NULLIFY (new_pos%rest) 625 CALL cp_sll_val_create(new_pos%rest, first_el=my_val) 626 new_pos => new_pos%rest 627 END IF 628 END IF 629 NULLIFY (my_val) 630 END IF 631 END DO 632 coord_section%values(ik, 1)%list => vals 633 634 CALL timestop(handle) 635 END SUBROUTINE section_coord_val_set 636 637! ************************************************************************************************** 638!> \brief routine to dump dipoles.. fast implementation 639!> \param dipoles ... 640!> \param dipoles_section ... 641!> \par History 642!> 12.2007 created [teo] 643!> \author Teodoro Laino 644! ************************************************************************************************** 645 SUBROUTINE update_dipoles_section(dipoles, dipoles_section) 646 647 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dipoles 648 TYPE(section_vals_type), POINTER :: dipoles_section 649 650 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_dipoles_section', & 651 routineP = moduleN//':'//routineN 652 653 INTEGER :: handle, ik, irk, Nlist, nloop 654 REAL(KIND=dp), DIMENSION(:), POINTER :: work 655 TYPE(cp_sll_val_type), POINTER :: new_pos, vals 656 TYPE(section_type), POINTER :: section 657 TYPE(val_type), POINTER :: my_val, old_val 658 659 CALL timeset(routineN, handle) 660 NULLIFY (my_val, old_val, section, vals) 661 CPASSERT(ASSOCIATED(dipoles_section)) 662 CPASSERT(dipoles_section%ref_count > 0) 663 section => dipoles_section%section 664 ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_") 665 IF (ik == -2) & 666 CALL cp_abort(__LOCATION__, & 667 "section "//TRIM(section%name)//" does not contain keyword "// & 668 "_DEFAULT_KEYWORD_") 669 670 ! At least one of the two arguments must be present.. 671 nloop = SIZE(dipoles, 2) 672 DO 673 IF (SIZE(dipoles_section%values, 2) == 1) EXIT 674 CALL section_vals_add_values(dipoles_section) 675 END DO 676 vals => dipoles_section%values(ik, 1)%list 677 Nlist = 0 678 IF (ASSOCIATED(vals)) THEN 679 Nlist = cp_sll_val_get_length(vals) 680 END IF 681 DO irk = 1, nloop 682 ALLOCATE (work(3)) 683 ! Always stored in A.U. 684 work = dipoles(1:3, irk) 685 CALL val_create(my_val, r_vals_ptr=work) 686 687 IF (Nlist /= 0) THEN 688 IF (irk == 1) THEN 689 new_pos => vals 690 ELSE 691 new_pos => new_pos%rest 692 END IF 693 old_val => new_pos%first_el 694 CALL val_release(old_val) 695 new_pos%first_el => my_val 696 ELSE 697 IF (irk == 1) THEN 698 NULLIFY (new_pos) 699 CALL cp_sll_val_create(new_pos, first_el=my_val) 700 vals => new_pos 701 ELSE 702 NULLIFY (new_pos%rest) 703 CALL cp_sll_val_create(new_pos%rest, first_el=my_val) 704 new_pos => new_pos%rest 705 END IF 706 END IF 707 NULLIFY (my_val) 708 END DO 709 dipoles_section%values(ik, 1)%list => vals 710 711 CALL timestop(handle) 712 END SUBROUTINE update_dipoles_section 713 714! ************************************************************************************************** 715!> \brief routine to dump quadrupoles.. fast implementation 716!> \param quadrupoles ... 717!> \param quadrupoles_section ... 718!> \par History 719!> 12.2007 created [teo] 720!> \author Teodoro Laino 721! ************************************************************************************************** 722 SUBROUTINE update_quadrupoles_section(quadrupoles, quadrupoles_section) 723 724 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: quadrupoles 725 TYPE(section_vals_type), POINTER :: quadrupoles_section 726 727 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_quadrupoles_section', & 728 routineP = moduleN//':'//routineN 729 730 INTEGER :: handle, i, ik, ind, irk, j, Nlist, nloop 731 REAL(KIND=dp), DIMENSION(:), POINTER :: work 732 TYPE(cp_sll_val_type), POINTER :: new_pos, vals 733 TYPE(section_type), POINTER :: section 734 TYPE(val_type), POINTER :: my_val, old_val 735 736 CALL timeset(routineN, handle) 737 NULLIFY (my_val, old_val, section, vals) 738 CPASSERT(ASSOCIATED(quadrupoles_section)) 739 CPASSERT(quadrupoles_section%ref_count > 0) 740 section => quadrupoles_section%section 741 ik = section_get_keyword_index(section, "_DEFAULT_KEYWORD_") 742 IF (ik == -2) & 743 CALL cp_abort(__LOCATION__, & 744 "section "//TRIM(section%name)//" does not contain keyword "// & 745 "_DEFAULT_KEYWORD_") 746 747 ! At least one of the two arguments must be present.. 748 nloop = SIZE(quadrupoles, 2) 749 DO 750 IF (SIZE(quadrupoles_section%values, 2) == 1) EXIT 751 CALL section_vals_add_values(quadrupoles_section) 752 END DO 753 vals => quadrupoles_section%values(ik, 1)%list 754 Nlist = 0 755 IF (ASSOCIATED(vals)) THEN 756 Nlist = cp_sll_val_get_length(vals) 757 END IF 758 DO irk = 1, nloop 759 ALLOCATE (work(6)) 760 ! Always stored in A.U. 761 ind = 0 762 DO i = 1, 3 763 DO j = i, 3 764 ind = ind + 1 765 work(ind) = quadrupoles(j, i, irk) 766 END DO 767 END DO 768 CALL val_create(my_val, r_vals_ptr=work) 769 770 IF (Nlist /= 0) THEN 771 IF (irk == 1) THEN 772 new_pos => vals 773 ELSE 774 new_pos => new_pos%rest 775 END IF 776 old_val => new_pos%first_el 777 CALL val_release(old_val) 778 new_pos%first_el => my_val 779 ELSE 780 IF (irk == 1) THEN 781 NULLIFY (new_pos) 782 CALL cp_sll_val_create(new_pos, first_el=my_val) 783 vals => new_pos 784 ELSE 785 NULLIFY (new_pos%rest) 786 CALL cp_sll_val_create(new_pos%rest, first_el=my_val) 787 new_pos => new_pos%rest 788 END IF 789 END IF 790 NULLIFY (my_val) 791 END DO 792 quadrupoles_section%values(ik, 1)%list => vals 793 794 CALL timestop(handle) 795 END SUBROUTINE update_quadrupoles_section 796 797! ************************************************************************************************** 798!> \brief Dump atomic, core, or shell coordinates to a file in CP2K &COORD 799!> section format 800!> \param particles ... 801!> \param molecules ... 802!> \param cell ... 803!> \param conv_factor ... 804!> \param output_unit ... 805!> \param core_or_shell ... 806!> \param scaled_coordinates ... 807!> \par History 808!> 07.02.2011 (Creation, MK) 809!> \author Matthias Krack (MK) 810!> \version 1.0 811! ************************************************************************************************** 812 SUBROUTINE dump_coordinates_cp2k(particles, molecules, cell, conv_factor, & 813 output_unit, core_or_shell, & 814 scaled_coordinates) 815 816 TYPE(particle_list_type), POINTER :: particles 817 TYPE(molecule_list_type), POINTER :: molecules 818 TYPE(cell_type), POINTER :: cell 819 REAL(KIND=dp), INTENT(IN) :: conv_factor 820 INTEGER, INTENT(IN) :: output_unit 821 LOGICAL, INTENT(IN) :: core_or_shell, scaled_coordinates 822 823 CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_coordinates_cp2k', & 824 routineP = moduleN//':'//routineN 825 826 CHARACTER(LEN=default_string_length) :: kind_name, tag 827 INTEGER :: handle, imolecule, iparticle, last_atom 828 LOGICAL :: ldum, molname_generated 829 REAL(KIND=dp), DIMENSION(3) :: r, s 830 TYPE(molecule_type), POINTER :: molecule 831 832 CALL timeset(routineN, handle) 833 834 CPASSERT(ASSOCIATED(particles)) 835 IF (.NOT. core_or_shell) THEN 836 CPASSERT(ASSOCIATED(molecules)) 837 END IF 838 IF (scaled_coordinates) THEN 839 CPASSERT(ASSOCIATED(cell)) 840 END IF 841 842 kind_name = "" 843 tag = "" 844 imolecule = 0 845 last_atom = 0 846 DO iparticle = 1, particles%n_els 847 CALL get_atomic_kind(particles%els(iparticle)%atomic_kind, name=kind_name) 848 IF (.NOT. core_or_shell) THEN 849 IF (iparticle > last_atom) THEN 850 imolecule = imolecule + 1 851 molecule => molecules%els(imolecule) 852 CALL get_molecule(molecule, last_atom=last_atom) 853 CALL get_molecule_kind(molecule%molecule_kind, & 854 molname_generated=molname_generated, & 855 name=tag) 856 IF (molname_generated) tag = "" 857 END IF 858 ldum = qmmm_ff_precond_only_qm(tag) 859 ldum = qmmm_ff_precond_only_qm(kind_name) 860 END IF 861 IF (scaled_coordinates) THEN 862 CALL real_to_scaled(s, particles%els(iparticle)%r, cell) 863 r(1:3) = s(1:3) 864 ELSE 865 r(1:3) = particles%els(iparticle)%r(1:3)*conv_factor 866 END IF 867 IF (core_or_shell) THEN 868 WRITE (UNIT=output_unit, FMT="(A,3(1X,ES25.16),1X,I0)") & 869 TRIM(ADJUSTL(kind_name)), r(1:3), particles%els(iparticle)%atom_index 870 ELSE 871 IF (LEN_TRIM(tag) > 0) THEN 872 WRITE (UNIT=output_unit, FMT="(A,3(1X,ES25.16),1X,A,1X,I0)") & 873 TRIM(ADJUSTL(kind_name)), r(1:3), TRIM(tag), imolecule 874 ELSE 875 WRITE (UNIT=output_unit, FMT="(A,3(1X,ES25.16))") & 876 TRIM(ADJUSTL(kind_name)), r(1:3) 877 END IF 878 END IF 879 END DO 880 881 CALL timestop(handle) 882 883 END SUBROUTINE dump_coordinates_cp2k 884 885END MODULE input_restart_force_eval 886