1!! Copyright (C) 2008 X. Andrade 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module ion_dynamics_oct_m 22 use iso_c_binding 23 use geometry_oct_m 24 use global_oct_m 25 use loct_math_oct_m 26 use messages_oct_m 27 use mpi_oct_m 28 use namespace_oct_m 29 use parser_oct_m 30 use read_coords_oct_m 31 use simul_box_oct_m 32 use species_oct_m 33 use tdfunction_oct_m 34 use profiling_oct_m 35 use unit_oct_m 36 use unit_system_oct_m 37 use varinfo_oct_m 38 39 implicit none 40 41 private 42 43 public :: & 44 ion_dynamics_t, & 45 ion_state_t, & 46 ion_dynamics_init, & 47 ion_dynamics_end, & 48 ion_dynamics_propagate, & 49 ion_dynamics_propagate_vel, & 50 ion_dynamics_save_state, & 51 ion_dynamics_restore_state, & 52 ion_dynamics_ions_move, & 53 ion_dynamics_temperature, & 54 ion_dynamics_kinetic_energy, & 55 ion_dynamics_freeze, & 56 ion_dynamics_unfreeze, & 57 ion_dynamics_verlet_step1, & 58 ion_dynamics_verlet_step2 59 60 integer, parameter :: & 61 THERMO_NONE = 0, & 62 THERMO_SCAL = 1, & 63 THERMO_NH = 2 64 65 type nose_hoover_t 66 private 67 FLOAT :: mass 68 FLOAT :: pos 69 FLOAT :: vel 70 end type nose_hoover_t 71 72 type ion_td_displacement_t 73 private 74 logical :: move 75 type(tdf_t) :: fx 76 type(tdf_t) :: fy 77 type(tdf_t) :: fz 78 end type ion_td_displacement_t 79 80 type ion_dynamics_t 81 private 82 logical :: move_ions 83 logical :: constant_velocity 84 integer :: thermostat 85 FLOAT :: dt 86 FLOAT :: current_temperature 87 88 FLOAT, pointer :: oldforce(:, :) 89 90 !> the old positions for Verlet (used for the Nose-Hoover) 91 FLOAT, pointer :: old_pos(:, :) 92 93 !> variables for the Nose-Hoover thermostat 94 type(nose_hoover_t) :: nh(1:2) 95 type(tdf_t) :: temperature_function 96 97 logical :: drive_ions 98 type(ion_td_displacement_t), pointer :: td_displacements(:) !> Time-dependent displacements driving the ions 99 type(geometry_t), pointer :: geo_t0 100 end type ion_dynamics_t 101 102 type ion_state_t 103 private 104 FLOAT, pointer :: pos(:, :) 105 FLOAT, pointer :: vel(:, :) 106 FLOAT, pointer :: old_pos(:, :) 107 type(nose_hoover_t) :: nh(1:2) 108 end type ion_state_t 109 110contains 111 112 ! --------------------------------------------------------- 113 subroutine ion_dynamics_init(this, namespace, geo) 114 type(ion_dynamics_t), intent(out) :: this 115 type(namespace_t), intent(in) :: namespace 116 type(geometry_t), intent(inout) :: geo 117 118 integer :: i, j, iatom, ierr 119 FLOAT :: x(MAX_DIM), temperature, sigma, kin1, kin2 120 type(c_ptr) :: random_gen_pointer 121 type(read_coords_info) :: xyz 122 character(len=100) :: temp_function_name 123 logical :: have_velocities 124 125 type(block_t) :: blk 126 integer :: ndisp 127 character(len=200) :: expression 128 129 PUSH_SUB(ion_dynamics_init) 130 131 nullify(this%oldforce) 132 133 have_velocities = .false. 134 this%drive_ions = .false. 135 136 !%Variable IonsConstantVelocity 137 !%Type logical 138 !%Default no 139 !%Section Time-Dependent::Propagation 140 !%Description 141 !% (Experimental) If this variable is set to yes, the ions will 142 !% move with a constant velocity given by the initial 143 !% conditions. They will not be affected by any forces. 144 !%End 145 call parse_variable(namespace, 'IonsConstantVelocity', .false., this%constant_velocity) 146 call messages_print_var_value(stdout, 'IonsConstantVelocity', this%constant_velocity) 147 148 if(this%constant_velocity) then 149 call messages_experimental('IonsConstantVelocity') 150 have_velocities = .true. 151 this%drive_ions = .true. 152 end if 153 154 !%Variable IonsTimeDependentDisplacements 155 !%Type block 156 !%Section Time-Dependent::Propagation 157 !%Description 158 !% (Experimental) This variable allows you to specify a 159 !% time-dependent function describing the displacement of the ions 160 !% from their equilibrium position: <math>r(t) = r_0 + \Delta 161 !% r(t)</math>. Specify the displacements dx(t), dy(t), dz(t) as 162 !% follows, for some or all of the atoms: 163 !% 164 !% <tt>%IonsTimeDependentDisplacements 165 !% <br> atom_index | "dx(t)" | "dy(t)" | "dz(t)" 166 !% <br>%</tt> 167 !% 168 !% The displacement functions are time-dependent functions and should match one 169 !% of the function names given in the first column of the <tt>TDFunctions</tt> block. 170 !% If this block is set, the ions will not be affected by any forces. 171 !%End 172 173 174 ndisp = 0 175 if(parse_block(namespace, 'IonsTimeDependentDisplacements', blk) == 0) then 176 call messages_experimental("IonsTimeDependentDisplacements") 177 ndisp= parse_block_n(blk) 178 SAFE_ALLOCATE(this%td_displacements(1:geo%natoms)) 179 this%td_displacements(1:geo%natoms)%move = .false. 180 if (ndisp > 0) this%drive_ions =.true. 181 182 do i = 1, ndisp 183 call parse_block_integer(blk, i-1, 0, iatom) 184 this%td_displacements(iatom)%move = .true. 185 186 call parse_block_string(blk, i-1, 1, expression) 187 call tdf_read(this%td_displacements(iatom)%fx, namespace, trim(expression), ierr) 188 if (ierr /= 0) then 189 write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:' 190 call messages_warning(1, namespace=namespace) 191 end if 192 193 194 call parse_block_string(blk, i-1, 2, expression) 195 call tdf_read(this%td_displacements(iatom)%fy, namespace, trim(expression), ierr) 196 if (ierr /= 0) then 197 write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:' 198 call messages_warning(1, namespace=namespace) 199 end if 200 201 call parse_block_string(blk, i-1, 3, expression) 202 call tdf_read(this%td_displacements(iatom)%fz, namespace, trim(expression), ierr) 203 if (ierr /= 0) then 204 write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:' 205 call messages_warning(1, namespace=namespace) 206 end if 207 208 end do 209 210 SAFE_ALLOCATE(this%geo_t0) 211 call geometry_copy(this%geo_t0, geo) 212 213 end if 214 215 216 217 218 !%Variable Thermostat 219 !%Type integer 220 !%Default none 221 !%Section Time-Dependent::Propagation 222 !%Description 223 !% This variable selects the type of thermostat applied to 224 !% control the ionic temperature. 225 !%Option none 0 226 !% No thermostat is applied. This is the default. 227 !%Option velocity_scaling 1 228 !% Velocities are scaled to control the temperature. 229 !%Option nose_hoover 2 230 !% Nose-Hoover thermostat. 231 !%End 232 233 call parse_variable(namespace, 'Thermostat', THERMO_NONE, this%thermostat) 234 if(.not.varinfo_valid_option('Thermostat', this%thermostat)) call messages_input_error(namespace, 'Thermostat') 235 call messages_print_var_option(stdout, 'Thermostat', this%thermostat) 236 237 if(this%thermostat /= THERMO_NONE) then 238 239 have_velocities = .true. 240 241 if(this%drive_ions) then 242 call messages_write('You cannot use a Thermostat and IonsConstantVelocity or IonsTimeDependentDisplacements') 243 call messages_write('at the same time.') 244 call messages_fatal(namespace=namespace) 245 end if 246 247 call messages_experimental('Thermostat') 248 249 !%Variable TemperatureFunction 250 !%Type integer 251 !%Default "temperature" 252 !%Section Time-Dependent::Propagation 253 !%Description 254 !% If a thermostat is used, this variable indicates the name of the 255 !% function in the <tt>TDFunctions</tt> block that will be used to control the 256 !% temperature. The values of the temperature are given in 257 !% degrees Kelvin. 258 !%End 259 call parse_variable(namespace, 'TemperatureFunction', 'temperature', temp_function_name) 260 261 call tdf_read(this%temperature_function, namespace, temp_function_name, ierr) 262 263 if(ierr /= 0) then 264 message(1) = "You have enabled a thermostat but Octopus could not find" 265 message(2) = "the '"//trim(temp_function_name)//"' function in the TDFunctions block." 266 call messages_fatal(2, namespace=namespace) 267 end if 268 269 if(this%thermostat == THERMO_NH) then 270 !%Variable ThermostatMass 271 !%Type float 272 !%Default 1.0 273 !%Section Time-Dependent::Propagation 274 !%Description 275 !% This variable sets the fictitious mass for the Nose-Hoover 276 !% thermostat. 277 !%End 278 call messages_obsolete_variable(namespace, 'NHMass', 'ThermostatMass') 279 280 call parse_variable(namespace, 'ThermostatMass', CNST(1.0), this%nh(1)%mass) 281 this%nh(2)%mass = this%nh(1)%mass 282 283 this%nh(1:2)%pos = M_ZERO 284 this%nh(1:2)%vel = M_ZERO 285 286 SAFE_ALLOCATE(this%old_pos(1:geo%space%dim, 1:geo%natoms)) 287 288 do iatom = 1, geo%natoms 289 this%old_pos(1:geo%space%dim, iatom) = geo%atom(iatom)%x(1:geo%space%dim) 290 end do 291 end if 292 293 end if 294 295 !now initialize velocities 296 297 !%Variable RandomVelocityTemp 298 !%Type float 299 !%Default 0.0 300 !%Section System::Velocities 301 !%Description 302 !% If this variable is present, <tt>Octopus</tt> will assign random 303 !% velocities to the atoms following a Boltzmann distribution with 304 !% temperature given by <tt>RandomVelocityTemp</tt> (in degrees Kelvin). 305 !%End 306 307 ! we now load the velocities, either from the temperature, from the input, or from a file 308 if(parse_is_defined(namespace, 'RandomVelocityTemp')) then 309 310 have_velocities = .true. 311 312 if( mpi_grp_is_root(mpi_world)) then 313 call loct_ran_init(random_gen_pointer) 314 call parse_variable(namespace, 'RandomVelocityTemp', M_ZERO, temperature, unit = unit_kelvin) 315 end if 316 317 do i = 1, geo%natoms 318 !generate the velocities in the root node 319 if( mpi_grp_is_root(mpi_world)) then 320 sigma = sqrt(temperature / species_mass(geo%atom(i)%species) ) 321 do j = 1, 3 322 geo%atom(i)%v(j) = loct_ran_gaussian(random_gen_pointer, sigma) 323 end do 324 end if 325#ifdef HAVE_MPI 326 !and send them to the others 327 call MPI_Bcast(geo%atom(i)%v, geo%space%dim, MPI_FLOAT, 0, mpi_world%comm, mpi_err) 328#endif 329 end do 330 331 if( mpi_grp_is_root(mpi_world)) then 332 call loct_ran_end(random_gen_pointer) 333 end if 334 335 kin1 = ion_dynamics_kinetic_energy(geo) 336 337 call cm_vel(geo, x) 338 do i = 1, geo%natoms 339 geo%atom(i)%v = geo%atom(i)%v - x 340 end do 341 342 kin2 = ion_dynamics_kinetic_energy(geo) 343 344 do i = 1, geo%natoms 345 geo%atom(i)%v(:) = sqrt(kin1/kin2)*geo%atom(i)%v(:) 346 end do 347 348 write(message(1),'(a,f10.4,1x,a)') 'Info: Initial velocities randomly distributed with T =', & 349 units_from_atomic(unit_kelvin, temperature), units_abbrev(unit_kelvin) 350 write(message(2),'(2x,a,f8.4,1x,a)') '<K> =', & 351 units_from_atomic(units_out%energy, ion_dynamics_kinetic_energy(geo)/geo%natoms), & 352 units_abbrev(units_out%energy) 353 write(message(3),'(2x,a,f8.4,1x,a)') '3/2 k_B T =', & 354 units_from_atomic(units_out%energy, (M_THREE/M_TWO)*temperature), & 355 units_abbrev(units_out%energy) 356 call messages_info(3) 357 358 else 359 !%Variable XYZVelocities 360 !%Type string 361 !%Section System::Velocities 362 !%Description 363 !% <tt>Octopus</tt> will try to read the starting velocities of the atoms from the XYZ file 364 !% specified by the variable <tt>XYZVelocities</tt>. 365 !% Note that you do not need to specify initial velocities if you are not going 366 !% to perform ion dynamics; if you are going to allow the ions to move but the velocities 367 !% are not specified, they are considered to be null. 368 !% Note: It is important for the velocities to maintain the ordering 369 !% in which the atoms were defined in the coordinates specifications. 370 !%End 371 372 !%Variable XSFVelocities 373 !%Type string 374 !%Section System::Velocities 375 !%Description 376 !% Like <tt>XYZVelocities</tt> but in XCrySDen format, as in <tt>XSFCoordinates</tt>. 377 !%End 378 379 !%Variable PDBVelocities 380 !%Type string 381 !%Section System::Velocities 382 !%Description 383 !% Like <tt>XYZVelocities</tt> but in PDB format, as in <tt>PDBCoordinates</tt>. 384 !%End 385 386 !%Variable Velocities 387 !%Type block 388 !%Section System::Velocities 389 !%Description 390 !% If <tt>XYZVelocities</tt>, <tt>PDBVelocities</tt>, and <tt>XSFVelocities</tt> 391 !% are not present, <tt>Octopus</tt> will try to fetch the initial 392 !% atomic velocities from this block. If this block is not present, <tt>Octopus</tt> 393 !% will set the initial velocities to zero. The format of this block can be 394 !% illustrated by this example: 395 !% 396 !% <tt>%Velocities 397 !% <br> 'C' | -1.7 | 0.0 | 0.0 398 !% <br> 'O' | 1.7 | 0.0 | 0.0 399 !% <br>%</tt> 400 !% 401 !% It describes one carbon and one oxygen moving at the relative 402 !% velocity of 3.4 velocity units. 403 !% 404 !% Note: It is important for the velocities to maintain the ordering 405 !% in which the atoms were defined in the coordinates specifications. 406 !%End 407 408 call read_coords_init(xyz) 409 call read_coords_read('Velocities', xyz, geo%space, namespace) 410 if(xyz%source /= READ_COORDS_ERR) then 411 412 have_velocities = .true. 413 414 if(geo%natoms /= xyz%n) then 415 write(message(1), '(a,i4,a,i4)') 'I need exactly ', geo%natoms, ' velocities, but I found ', xyz%n 416 call messages_fatal(1, namespace=namespace) 417 end if 418 419 ! copy information and adjust units 420 do i = 1, geo%natoms 421 geo%atom(i)%v = units_to_atomic(units_inp%velocity/units_inp%length, xyz%atom(i)%x) 422 end do 423 424 call read_coords_end(xyz) 425 426 else 427 do i = 1, geo%natoms 428 geo%atom(i)%v = M_ZERO 429 end do 430 end if 431 end if 432 433 geo%kinetic_energy = ion_dynamics_kinetic_energy(geo) 434 435 !%Variable MoveIons 436 !%Type logical 437 !%Section Time-Dependent::Propagation 438 !%Description 439 !% This variable controls whether atoms are moved during a time 440 !% propagation run. The default is yes when the ion velocity is 441 !% set explicitly or implicitly, otherwise is no. 442 !%End 443 call parse_variable(namespace, 'MoveIons', have_velocities, this%move_ions) 444 call messages_print_var_value(stdout, 'MoveIons', this%move_ions) 445 446 if(ion_dynamics_ions_move(this)) then 447 SAFE_ALLOCATE(this%oldforce(1:geo%space%dim, 1:geo%natoms)) 448 end if 449 450 POP_SUB(ion_dynamics_init) 451 end subroutine ion_dynamics_init 452 453 454 ! --------------------------------------------------------- 455 subroutine ion_dynamics_end(this) 456 type(ion_dynamics_t), intent(inout) :: this 457 458 PUSH_SUB(ion_dynamics_end) 459 SAFE_DEALLOCATE_P(this%oldforce) 460 461 if(this%thermostat /= THERMO_NONE) then 462 call tdf_end(this%temperature_function) 463 end if 464 465 if (this%drive_ions .and. associated(this%td_displacements) ) then 466 if (any (this%td_displacements(1:this%geo_t0%natoms)%move)) then 467 ! geometry end cannot be called here, otherwise the species are destroyed twice 468 ! call geometry_end(this%geo_t0) 469 end if 470 SAFE_DEALLOCATE_P(this%td_displacements) 471 if (any (this%td_displacements(:)%move)) then 472 call geometry_end(this%geo_t0) 473 end if 474 end if 475 476 477 478 479 POP_SUB(ion_dynamics_end) 480 end subroutine ion_dynamics_end 481 482 483 ! --------------------------------------------------------- 484 subroutine ion_dynamics_propagate(this, sb, geo, time, dt, namespace) 485 type(ion_dynamics_t), intent(inout) :: this 486 type(simul_box_t), intent(in) :: sb 487 type(geometry_t), intent(inout) :: geo 488 FLOAT, intent(in) :: time 489 FLOAT, intent(in) :: dt 490 type(namespace_t), intent(in) :: namespace 491 492 integer :: iatom 493 FLOAT :: DR(1:3) 494 495 if(.not. ion_dynamics_ions_move(this)) return 496 497 PUSH_SUB(ion_dynamics_propagate) 498 499 DR = M_ZERO 500 501 this%dt = dt 502 503 504 ! get the temperature from the tdfunction for the current time 505 if(this%thermostat /= THERMO_NONE) then 506 this%current_temperature = units_to_atomic(unit_kelvin, tdf(this%temperature_function, time)) 507 508 if(this%current_temperature < M_ZERO) then 509 write(message(1), '(a, f10.3, 3a, f10.3, 3a)') & 510 "Negative temperature (", & 511 units_from_atomic(unit_kelvin, this%current_temperature), " ", units_abbrev(unit_kelvin), & 512 ") at time ", & 513 units_from_atomic(units_out%time, time), " ", trim(units_abbrev(units_out%time)), "." 514 call messages_fatal(1, namespace=namespace) 515 end if 516 else 517 this%current_temperature = CNST(0.0) 518 end if 519 520 if(this%thermostat /= THERMO_NH) then 521 ! integrate using verlet 522 do iatom = 1, geo%natoms 523 if(.not. geo%atom(iatom)%move) cycle 524 525 if(.not. this%drive_ions) then 526 527 geo%atom(iatom)%x(1:geo%space%dim) = geo%atom(iatom)%x(1:geo%space%dim) & 528 + dt*geo%atom(iatom)%v(1:geo%space%dim) + & 529 M_HALF*dt**2 / species_mass(geo%atom(iatom)%species) * geo%atom(iatom)%f(1:geo%space%dim) 530 531 this%oldforce(1:geo%space%dim, iatom) = geo%atom(iatom)%f(1:geo%space%dim) 532 533 else 534 if(this%constant_velocity) then 535 geo%atom(iatom)%x(1:geo%space%dim) = geo%atom(iatom)%x(1:geo%space%dim) & 536 + dt*geo%atom(iatom)%v(1:geo%space%dim) 537 end if 538 539 540 if (this%td_displacements(iatom)%move) then 541 542 DR(1:3)=(/TOFLOAT(tdf(this%td_displacements(iatom)%fx,time)), & 543 TOFLOAT(tdf(this%td_displacements(iatom)%fy,time)), & 544 TOFLOAT(tdf(this%td_displacements(iatom)%fz,time)) /) 545 546 geo%atom(iatom)%x(1:geo%space%dim) = this%geo_t0%atom(iatom)%x(1:geo%space%dim) + DR(1:geo%space%dim) 547 end if 548 549 end if 550 551 end do 552 553 else 554 ! for the Nose-Hoover thermostat we use a special integrator 555 556 ! The implementation of the Nose-Hoover thermostat is based on 557 ! Understanding Molecular Simulations by Frenkel and Smit, 558 ! Appendix E, page 540-542. 559 560 call nh_chain(this, geo) 561 562 do iatom = 1, geo%natoms 563 geo%atom(iatom)%x(1:geo%space%dim) = geo%atom(iatom)%x(1:geo%space%dim) + M_HALF*dt*geo%atom(iatom)%v(1:geo%space%dim) 564 end do 565 566 end if 567 568 call simul_box_atoms_in_box(sb, geo, namespace, .false.) 569 570 POP_SUB(ion_dynamics_propagate) 571 end subroutine ion_dynamics_propagate 572 573 574 ! --------------------------------------------------------- 575 subroutine nh_chain(this, geo) 576 type(ion_dynamics_t), intent(inout) :: this 577 type(geometry_t), intent(inout) :: geo 578 579 FLOAT :: g1, g2, ss, uk, dt, temp 580 integer :: iatom 581 582 PUSH_SUB(nh_chain) 583 584 dt = this%dt 585 586 uk = ion_dynamics_kinetic_energy(geo) 587 588 temp = this%current_temperature 589 590 g2 = (this%nh(1)%mass*this%nh(1)%vel**2 - temp)/this%nh(2)%mass 591 this%nh(2)%vel = this%nh(2)%vel + g2*dt/CNST(4.0) 592 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/CNST(8.0)) 593 594 g1 = (CNST(2.0)*uk - M_THREE*geo%natoms*temp)/this%nh(1)%mass 595 this%nh(1)%vel = this%nh(1)%vel + g1*dt/CNST(4.0) 596 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/CNST(8.0)) 597 this%nh(1)%pos = this%nh(1)%pos + this%nh(1)%vel*dt/CNST(2.0) 598 this%nh(2)%pos = this%nh(2)%pos + this%nh(2)%vel*dt/CNST(2.0) 599 600 ss = exp(-this%nh(1)%vel*dt/CNST(2.0)) 601 602 do iatom = 1, geo%natoms 603 geo%atom(iatom)%v(1:geo%space%dim) = ss*geo%atom(iatom)%v(1:geo%space%dim) 604 end do 605 606 uk = uk*ss**2 607 608 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/CNST(8.0)) 609 g1 = (CNST(2.0)*uk - M_THREE*geo%natoms*temp)/this%nh(1)%mass 610 this%nh(1)%vel = this%nh(1)%vel + g1*dt/CNST(4.0) 611 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/CNST(8.0)) 612 613 g2 = (this%nh(1)%mass*this%nh(1)%vel**2 - temp)/this%nh(2)%mass 614 this%nh(2)%vel = this%nh(2)%vel + g2*dt/CNST(4.0) 615 616 POP_SUB(nh_chain) 617 end subroutine nh_chain 618 619 620 ! --------------------------------------------------------- 621 subroutine ion_dynamics_propagate_vel(this, geo, atoms_moved) 622 type(ion_dynamics_t), intent(inout) :: this 623 type(geometry_t), intent(inout) :: geo 624 logical, optional, intent(out) :: atoms_moved !< Returns true if the atoms were moved by this function. 625 626 integer :: iatom 627 FLOAT :: scal 628 629 if(.not. ion_dynamics_ions_move(this)) return 630 if(this%drive_ions) return 631 632 PUSH_SUB(ion_dynamics_propagate_vel) 633 634 if(present(atoms_moved)) atoms_moved = this%thermostat == THERMO_NH 635 636 if(this%thermostat /= THERMO_NH) then 637 ! velocity verlet 638 639 do iatom = 1, geo%natoms 640 if(.not. geo%atom(iatom)%move) cycle 641 642 geo%atom(iatom)%v(1:geo%space%dim) = geo%atom(iatom)%v(1:geo%space%dim) & 643 + this%dt/species_mass(geo%atom(iatom)%species) * M_HALF * (this%oldforce(1:geo%space%dim, iatom) + & 644 geo%atom(iatom)%f(1:geo%space%dim)) 645 646 end do 647 648 else 649 ! the nose-hoover integration 650 do iatom = 1, geo%natoms 651 geo%atom(iatom)%v(1:geo%space%dim) = geo%atom(iatom)%v(1:geo%space%dim) + & 652 this%dt*geo%atom(iatom)%f(1:geo%space%dim) / species_mass(geo%atom(iatom)%species) 653 geo%atom(iatom)%x(1:geo%space%dim) = geo%atom(iatom)%x(1:geo%space%dim) + M_HALF*this%dt*geo%atom(iatom)%v(1:geo%space%dim) 654 end do 655 656 call nh_chain(this, geo) 657 658 end if 659 660 if(this%thermostat == THERMO_SCAL) then 661 scal = sqrt(this%current_temperature/ion_dynamics_temperature(geo)) 662 663 do iatom = 1, geo%natoms 664 geo%atom(iatom)%v(1:geo%space%dim) = scal*geo%atom(iatom)%v(1:geo%space%dim) 665 end do 666 end if 667 668 POP_SUB(ion_dynamics_propagate_vel) 669 end subroutine ion_dynamics_propagate_vel 670 671 672 ! --------------------------------------------------------- 673 !> A bare verlet integrator. 674 subroutine ion_dynamics_verlet_step1(geo, q, v, fold, dt) 675 type(geometry_t), intent(in) :: geo 676 FLOAT, intent(inout) :: q(:, :) 677 FLOAT, intent(inout) :: v(:, :) 678 FLOAT, intent(in) :: fold(:, :) 679 FLOAT, intent(in) :: dt 680 681 integer :: iatom 682 683 PUSH_SUB(ion_dynamics_verlet_step1) 684 685 ! First transform momenta to velocities 686 do iatom = 1, geo%natoms 687 v(iatom, 1:geo%space%dim) = v(iatom, 1:geo%space%dim) / species_mass(geo%atom(iatom)%species) 688 end do 689 690 ! integrate using verlet 691 do iatom = 1, geo%natoms 692 if(.not. geo%atom(iatom)%move) cycle 693 q(iatom, 1:geo%space%dim) = q(iatom, 1:geo%space%dim) & 694 + dt * v(iatom, 1:geo%space%dim) + & 695 M_HALF*dt**2 / species_mass(geo%atom(iatom)%species) * fold(iatom, 1:geo%space%dim) 696 end do 697 698 ! And back to momenta. 699 do iatom = 1, geo%natoms 700 v(iatom, 1:geo%space%dim) = species_mass(geo%atom(iatom)%species) * v(iatom, 1:geo%space%dim) 701 end do 702 703 POP_SUB(ion_dynamics_verlet_step1) 704 end subroutine ion_dynamics_verlet_step1 705 706 707 708 ! --------------------------------------------------------- 709 !> A bare verlet integrator. 710 subroutine ion_dynamics_verlet_step2(geo, v, fold, fnew, dt) 711 type(geometry_t), intent(in) :: geo 712 FLOAT, intent(inout) :: v(:, :) 713 FLOAT, intent(in) :: fold(:, :) 714 FLOAT, intent(in) :: fnew(:, :) 715 FLOAT, intent(in) :: dt 716 717 integer :: iatom 718 719 PUSH_SUB(ion_dynamics_verlet_step2) 720 721 ! First transform momenta to velocities 722 do iatom = 1, geo%natoms 723 v(iatom, 1:geo%space%dim) = v(iatom, 1:geo%space%dim) / species_mass(geo%atom(iatom)%species) 724 end do 725 726 ! velocity verlet 727 do iatom = 1, geo%natoms 728 if(.not. geo%atom(iatom)%move) cycle 729 v(iatom, 1:geo%space%dim) = v(iatom, 1:geo%space%dim) & 730 + dt / species_mass(geo%atom(iatom)%species) * M_HALF * (fold(iatom, 1:geo%space%dim) + & 731 fnew(iatom, 1:geo%space%dim)) 732 end do 733 734 ! And back to momenta. 735 do iatom = 1, geo%natoms 736 v(iatom, 1:geo%space%dim) = species_mass(geo%atom(iatom)%species) * v(iatom, 1:geo%space%dim) 737 end do 738 739 POP_SUB(ion_dynamics_verlet_step2) 740 end subroutine ion_dynamics_verlet_step2 741 742 743 ! --------------------------------------------------------- 744 subroutine ion_dynamics_save_state(this, geo, state) 745 type(ion_dynamics_t), intent(in) :: this 746 type(geometry_t), intent(in) :: geo 747 type(ion_state_t), intent(out) :: state 748 749 integer :: iatom 750 751 if(.not. ion_dynamics_ions_move(this)) return 752 753 PUSH_SUB(ion_dynamics_save_state) 754 755 SAFE_ALLOCATE(state%pos(1:geo%space%dim, 1:geo%natoms)) 756 SAFE_ALLOCATE(state%vel(1:geo%space%dim, 1:geo%natoms)) 757 758 do iatom = 1, geo%natoms 759 state%pos(1:geo%space%dim, iatom) = geo%atom(iatom)%x(1:geo%space%dim) 760 state%vel(1:geo%space%dim, iatom) = geo%atom(iatom)%v(1:geo%space%dim) 761 end do 762 763 if(this%thermostat == THERMO_NH) then 764 SAFE_ALLOCATE(state%old_pos(1:geo%space%dim, 1:geo%natoms)) 765 state%old_pos(1:geo%space%dim, 1:geo%natoms) = this%old_pos(1:geo%space%dim, 1:geo%natoms) 766 state%nh(1:2)%pos = this%nh(1:2)%pos 767 state%nh(1:2)%vel = this%nh(1:2)%vel 768 end if 769 770 POP_SUB(ion_dynamics_save_state) 771 end subroutine ion_dynamics_save_state 772 773 774 ! --------------------------------------------------------- 775 subroutine ion_dynamics_restore_state(this, geo, state) 776 type(ion_dynamics_t), intent(inout) :: this 777 type(geometry_t), intent(inout) :: geo 778 type(ion_state_t), intent(inout) :: state 779 780 integer :: iatom 781 782 if(.not. ion_dynamics_ions_move(this)) return 783 784 PUSH_SUB(ion_dynamics_restore_state) 785 786 do iatom = 1, geo%natoms 787 geo%atom(iatom)%x(1:geo%space%dim) = state%pos(1:geo%space%dim, iatom) 788 geo%atom(iatom)%v(1:geo%space%dim) = state%vel(1:geo%space%dim, iatom) 789 end do 790 791 if(this%thermostat == THERMO_NH) then 792 this%old_pos(1:geo%space%dim, 1:geo%natoms) = state%old_pos(1:geo%space%dim, 1:geo%natoms) 793 this%nh(1:2)%pos = state%nh(1:2)%pos 794 this%nh(1:2)%vel = state%nh(1:2)%vel 795 SAFE_DEALLOCATE_P(state%old_pos) 796 end if 797 798 SAFE_DEALLOCATE_P(state%pos) 799 SAFE_DEALLOCATE_P(state%vel) 800 801 POP_SUB(ion_dynamics_restore_state) 802 end subroutine ion_dynamics_restore_state 803 804 805 ! --------------------------------------------------------- 806 logical pure function ion_dynamics_ions_move(this) result(ions_move) 807 type(ion_dynamics_t), intent(in) :: this 808 809 ions_move = this%move_ions 810 811 end function ion_dynamics_ions_move 812 813 814 ! --------------------------------------------------------- 815 FLOAT pure function ion_dynamics_kinetic_energy(geo) result(kinetic_energy) 816 type(geometry_t), intent(in) :: geo 817 818 integer :: iatom 819 820 kinetic_energy = M_ZERO 821 do iatom = 1, geo%natoms 822 kinetic_energy = kinetic_energy + & 823 M_HALF * species_mass(geo%atom(iatom)%species) * sum(geo%atom(iatom)%v(1:geo%space%dim)**2) 824 end do 825 826 end function ion_dynamics_kinetic_energy 827 828 829 ! --------------------------------------------------------- 830 !> This function returns the ionic temperature in energy units. 831 FLOAT pure function ion_dynamics_temperature(geo) result(temperature) 832 type(geometry_t), intent(in) :: geo 833 834 temperature = CNST(2.0)/CNST(3.0)*ion_dynamics_kinetic_energy(geo)/geo%natoms 835 836 end function ion_dynamics_temperature 837 838 839 ! --------------------------------------------------------- 840 !> Freezes the ionic movement. 841 logical function ion_dynamics_freeze(this) result(freeze) 842 type(ion_dynamics_t), intent(inout) :: this 843 if(this%move_ions) then 844 this%move_ions = .false. 845 freeze = .true. 846 else 847 freeze = .false. 848 end if 849 end function ion_dynamics_freeze 850 851 852 ! --------------------------------------------------------- 853 !> Unfreezes the ionic movement. 854 subroutine ion_dynamics_unfreeze(this) 855 type(ion_dynamics_t), intent(inout) :: this 856 this%move_ions = .true. 857 end subroutine ion_dynamics_unfreeze 858 859 860end module ion_dynamics_oct_m 861 862!! Local Variables: 863!! mode: f90 864!! coding: utf-8 865!! End: 866