Octopus
ion_dynamics.F90
Go to the documentation of this file.
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
22 use accel_oct_m
24 use iso_c_binding
25 use debug_oct_m
26 use global_oct_m
27 use grid_oct_m
28 use ions_oct_m
29 use, intrinsic :: iso_fortran_env
34 use math_oct_m
36 use mpi_oct_m
39 use parser_oct_m
43 use space_oct_m
48 use unit_oct_m
51
52 implicit none
53
54 private
55
56 public :: &
74
75
76 integer, parameter :: &
77 THERMO_NONE = 0, &
78 thermo_scal = 1, &
79 thermo_nh = 2
80
81 type nose_hoover_t
82 private
83 real(real64) :: mass
84 real(real64) :: pos
85 real(real64) :: vel
86 end type nose_hoover_t
87
89 private
90 logical :: move
91 type(tdf_t) :: fx
92 type(tdf_t) :: fy
93 type(tdf_t) :: fz
95
97 private
98 logical :: move_ions
99 logical :: relax_cell
100 logical :: constant_velocity
101 integer :: thermostat
102 real(real64) :: dt
103 real(real64) :: current_temperature
104
105 real(real64), allocatable :: oldforce(:, :)
106
108 real(real64), allocatable :: old_pos(:, :)
109
111 real(real64), allocatable :: cell_force(:)
112 real(real64), allocatable :: old_cell_force(:)
113 real(real64), allocatable :: cell_vel(:)
114 real(real64), allocatable :: initial_rlattice(:,:)
115 real(real64), allocatable :: strain(:)
116
117 real(real64) :: pressure
118
120 logical :: symmetrize = .false.
121 type(symmetrizer_t), pointer :: symm
122
124 type(nose_hoover_t) :: nh(1:2)
125 type(tdf_t) :: temperature_function
126
128 logical :: drive_ions
129 type(ion_td_displacement_t), allocatable :: td_displacements(:)
130 type(ions_t), pointer :: ions_t0
131
132 real(real64), public :: ionic_scale
133 contains
134 procedure :: ions_move => ion_dynamics_ions_move
135 procedure :: cell_relax => ion_dynamics_cell_relax
136 procedure :: update_stress => ion_dynamics_update_stress
137 procedure :: is_active => ion_dynamics_is_active
138 end type ion_dynamics_t
139
140 type ion_state_t
141 private
142 real(real64), allocatable :: pos(:, :)
143 real(real64), allocatable :: vel(:, :)
144 real(real64), allocatable :: old_pos(:, :)
145 type(nose_hoover_t) :: nh(1:2)
146 end type ion_state_t
147
148 type cell_state_t
149 private
150 real(real64), allocatable :: pos(:, :)
151 real(real64), allocatable :: vel(:, :)
152 real(real64), allocatable :: old_pos(:, :)
153 end type cell_state_t
154
155contains
156
157 ! ---------------------------------------------------------
158 subroutine ion_dynamics_init(this, namespace, ions, symmetrize, symm)
159 use iso_c_binding, only: c_ptr
160 type(ion_dynamics_t), intent(out) :: this
161 type(namespace_t), intent(in) :: namespace
162 type(ions_t), intent(inout) :: ions
163 logical, intent(in) :: symmetrize
164 type(symmetrizer_t), optional, target, intent(in) :: symm
165
166 integer :: i, j, iatom, ierr, periodic_dim, ncomp
167 real(real64) :: xx(ions%space%dim), temperature, sigma, kin1, kin2
168 type(c_ptr) :: random_gen_pointer
169 type(read_coords_info) :: xyz
170 character(len=100) :: temp_function_name
171 logical :: have_velocities
172
173 type(block_t) :: blk
174 integer :: ndisp
175 character(len=200) :: expression
179 have_velocities = .false.
180 this%drive_ions = .false.
182 this%symmetrize = symmetrize
183 if (this%symmetrize) then
184 assert(present(symm))
185 this%symm => symm
186 else
187 nullify(this%symm)
188 end if
190 !%Variable TDIonicTimeScale
191 !%Type float
192 !%Default 1.0
193 !%Section Time-Dependent::Propagation
194 !%Description
195 !% This variable defines the factor between the timescale of ionic
196 !% and electronic movement. It allows reasonably fast
197 !% Born-Oppenheimer molecular-dynamics simulations based on
198 !% Ehrenfest dynamics. The value of this variable is equivalent to
199 !% the role of <math>\mu</math> in Car-Parrinello. Increasing it
200 !% linearly accelerates the time step of the ion
201 !% dynamics, but also increases the deviation of the system from the
202 !% Born-Oppenheimer surface. The default is 1, which means that both
203 !% timescales are the same. Note that a value different than 1
204 !% implies that the electrons will not follow physical behaviour.
205 !%
206 !% According to our tests, values around 10 are reasonable, but it
207 !% will depend on your system, mainly on the width of the gap.
208 !%
209 !% Important: The electronic time step will be the value of
210 !% <tt>TDTimeStep</tt> divided by this variable, so if you have determined an
211 !% optimal electronic time step (that we can call <i>dte</i>), it is
212 !% recommended that you define your time step as:
213 !%
214 !% <tt>TDTimeStep</tt> = <i>dte</i> * <tt>TDIonicTimeScale</tt>
215 !%
216 !% so you will always use the optimal electronic time step
217 !% (<a href=http://arxiv.org/abs/0710.3321>more details</a>).
218 !%End
219 call parse_variable(namespace, 'TDIonicTimeScale', m_one, this%ionic_scale)
220
221 if (this%ionic_scale <= m_zero) then
222 write(message(1),'(a)') 'Input: TDIonicTimeScale must be positive.'
223 call messages_fatal(1, namespace=namespace)
224 end if
226 call messages_print_var_value('TDIonicTimeScale', this%ionic_scale, namespace=namespace)
231 !%Variable IonsConstantVelocity
232 !%Type logical
233 !%Default no
234 !%Section Time-Dependent::Propagation
235 !%Description
236 !% (Experimental) If this variable is set to yes, the ions will
237 !% move with a constant velocity given by the initial
238 !% conditions. They will not be affected by any forces.
239 !%End
240 call parse_variable(namespace, 'IonsConstantVelocity', .false., this%constant_velocity)
241 call messages_print_var_value('IonsConstantVelocity', this%constant_velocity, namespace=namespace)
242
243 if (this%constant_velocity) then
244 call messages_experimental('IonsConstantVelocity', namespace=namespace)
245 have_velocities = .true.
246 this%drive_ions = .true.
247 end if
248
249 !%Variable IonsTimeDependentDisplacements
250 !%Type block
251 !%Section Time-Dependent::Propagation
252 !%Description
253 !% (Experimental) This variable allows you to specify a
254 !% time-dependent function describing the displacement of the ions
255 !% from their equilibrium position: <math>r(t) = r_0 + \Delta
256 !% r(t)</math>. Specify the displacements dx(t), dy(t), dz(t) as
257 !% follows, for some or all of the atoms:
258 !%
259 !% <tt>%IonsTimeDependentDisplacements
260 !% <br>&nbsp;&nbsp; atom_index | "dx(t)" | "dy(t)" | "dz(t)"
261 !% <br>%</tt>
262 !%
263 !% The displacement functions are time-dependent functions and should match one
264 !% of the function names given in the first column of the <tt>TDFunctions</tt> block.
265 !% If this block is set, the ions will not be affected by any forces.
266 !%End
267
268
269 ndisp = 0
270 if (parse_block(namespace, 'IonsTimeDependentDisplacements', blk) == 0) then
271 call messages_experimental("IonsTimeDependentDisplacements", namespace=namespace)
272 ndisp= parse_block_n(blk)
273 safe_allocate(this%td_displacements(1:ions%natoms))
274 this%td_displacements(1:ions%natoms)%move = .false.
275 if (ndisp > 0) this%drive_ions =.true.
276
277 do i = 1, ndisp
278 call parse_block_integer(blk, i-1, 0, iatom)
279 this%td_displacements(iatom)%move = .true.
280
281 call parse_block_string(blk, i-1, 1, expression)
282 call tdf_read(this%td_displacements(iatom)%fx, namespace, trim(expression), ierr)
283 if (ierr /= 0) then
284 write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:'
285 call messages_warning(1, namespace=namespace)
286 end if
287
288
289 call parse_block_string(blk, i-1, 2, expression)
290 call tdf_read(this%td_displacements(iatom)%fy, namespace, trim(expression), ierr)
291 if (ierr /= 0) then
292 write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:'
293 call messages_warning(1, namespace=namespace)
294 end if
295
296 call parse_block_string(blk, i-1, 3, expression)
297 call tdf_read(this%td_displacements(iatom)%fz, namespace, trim(expression), ierr)
298 if (ierr /= 0) then
299 write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:'
300 call messages_warning(1, namespace=namespace)
301 end if
302
303 end do
304
305 safe_allocate(this%ions_t0)
306 this%ions_t0 = ions
307
308 end if
309
310
311
312
313 !%Variable Thermostat
314 !%Type integer
315 !%Default none
316 !%Section Time-Dependent::Propagation
317 !%Description
318 !% This variable selects the type of thermostat applied to
319 !% control the ionic temperature.
320 !%Option none 0
321 !% No thermostat is applied. This is the default.
322 !%Option velocity_scaling 1
323 !% Velocities are scaled to control the temperature.
324 !%Option nose_hoover 2
325 !% Nose-Hoover thermostat.
326 !%End
327
328 call parse_variable(namespace, 'Thermostat', thermo_none, this%thermostat)
329 if (.not. varinfo_valid_option('Thermostat', this%thermostat)) call messages_input_error(namespace, 'Thermostat')
330 call messages_print_var_option('Thermostat', this%thermostat, namespace=namespace)
331
332 if (this%thermostat /= thermo_none) then
333
334 have_velocities = .true.
335
336 if (this%drive_ions) then
337 call messages_write('You cannot use a Thermostat and IonsConstantVelocity or IonsTimeDependentDisplacements')
338 call messages_write('at the same time.')
339 call messages_fatal(namespace=namespace)
340 end if
341
342 call messages_experimental('Thermostat', namespace=namespace)
343
344 !%Variable TemperatureFunction
345 !%Type integer
346 !%Default "temperature"
347 !%Section Time-Dependent::Propagation
348 !%Description
349 !% If a thermostat is used, this variable indicates the name of the
350 !% function in the <tt>TDFunctions</tt> block that will be used to control the
351 !% temperature. The values of the temperature are given in
352 !% degrees Kelvin.
353 !%End
354 call parse_variable(namespace, 'TemperatureFunction', 'temperature', temp_function_name)
355
356 call tdf_read(this%temperature_function, namespace, temp_function_name, ierr)
357
358 if (ierr /= 0) then
359 message(1) = "You have enabled a thermostat but Octopus could not find"
360 message(2) = "the '"//trim(temp_function_name)//"' function in the TDFunctions block."
361 call messages_fatal(2, namespace=namespace)
362 end if
363
364 if (this%thermostat == thermo_nh) then
365 !%Variable ThermostatMass
366 !%Type float
367 !%Default 1.0
368 !%Section Time-Dependent::Propagation
369 !%Description
370 !% This variable sets the fictitious mass for the Nose-Hoover
371 !% thermostat.
372 !%End
373 call messages_obsolete_variable(namespace, 'NHMass', 'ThermostatMass')
374
375 call parse_variable(namespace, 'ThermostatMass', m_one, this%nh(1)%mass)
376 this%nh(2)%mass = this%nh(1)%mass
377
378 this%nh(1:2)%pos = m_zero
379 this%nh(1:2)%vel = m_zero
380
381 safe_allocate(this%old_pos(1:ions%space%dim, 1:ions%natoms))
382
383 this%old_pos = ions%pos
384 end if
385
386 end if
387
388 !now initialize velocities
389
390 !%Variable RandomVelocityTemp
391 !%Type float
392 !%Default 0.0
393 !%Section System::Velocities
394 !%Description
395 !% If this variable is present, <tt>Octopus</tt> will assign random
396 !% velocities to the atoms following a Boltzmann distribution with
397 !% temperature given by <tt>RandomVelocityTemp</tt> (in degrees Kelvin).
398 !% The seed for the random number generator can be modified by setting
399 !% <tt>GSL_RNG_SEED</tt> environment variable.
400 !%End
401
402 ! we now load the velocities, either from the temperature, from the input, or from a file
403 if (parse_is_defined(namespace, 'RandomVelocityTemp')) then
404
405 have_velocities = .true.
406
407 if (mpi_grp_is_root(mpi_world)) then
408 call loct_ran_init(random_gen_pointer)
409 call parse_variable(namespace, 'RandomVelocityTemp', m_zero, temperature, unit = unit_kelvin)
410 end if
411
412 do i = 1, ions%natoms
413 !generate the velocities in the root node
414 if (mpi_grp_is_root(mpi_world)) then
415 sigma = sqrt(temperature / ions%mass(i))
416 do j = 1, 3
417 ions%vel(j, i) = loct_ran_gaussian(random_gen_pointer, sigma)
418 end do
419 end if
420 !and send them to the others
421 call mpi_world%bcast(ions%vel(:, i), ions%space%dim, mpi_double_precision, 0)
422 end do
423
424 if (mpi_grp_is_root(mpi_world)) then
425 call loct_ran_end(random_gen_pointer)
426 end if
427
428 call ions%update_kinetic_energy()
429 kin1 = ions%kinetic_energy
430
431 xx = ions%center_of_mass_vel()
432 do i = 1, ions%natoms
433 ions%vel(:, i) = ions%vel(:, i) - xx
434 end do
435
436 call ions%update_kinetic_energy()
437 kin2 = ions%kinetic_energy
438
439 do i = 1, ions%natoms
440 ions%vel(:, i) = sqrt(kin1/kin2)*ions%vel(:, i)
441 end do
442
443 call ions%update_kinetic_energy()
444
445 write(message(1),'(a,f10.4,1x,a)') 'Info: Initial velocities randomly distributed with T =', &
447 write(message(2),'(2x,a,f8.4,1x,a)') '<K> =', &
448 units_from_atomic(units_out%energy, ions%kinetic_energy/ions%natoms), &
449 units_abbrev(units_out%energy)
450 write(message(3),'(2x,a,f8.4,1x,a)') '3/2 k_B T =', &
451 units_from_atomic(units_out%energy, (m_three/m_two)*temperature), &
452 units_abbrev(units_out%energy)
453 call messages_info(3, namespace=namespace)
454
455 else
456 !%Variable XYZVelocities
457 !%Type string
458 !%Section System::Velocities
459 !%Description
460 !% <tt>Octopus</tt> will try to read the starting velocities of the atoms from the XYZ file
461 !% specified by the variable <tt>XYZVelocities</tt>.
462 !% Note that you do not need to specify initial velocities if you are not going
463 !% to perform ion dynamics; if you are going to allow the ions to move but the velocities
464 !% are not specified, they are considered to be null.
465 !% Note: It is important for the velocities to maintain the ordering
466 !% in which the atoms were defined in the coordinates specifications.
467 !%End
468
469 !%Variable XSFVelocities
470 !%Type string
471 !%Section System::Velocities
472 !%Description
473 !% Like <tt>XYZVelocities</tt> but in XCrySDen format, as in <tt>XSFCoordinates</tt>.
474 !%End
475
476 !%Variable PDBVelocities
477 !%Type string
478 !%Section System::Velocities
479 !%Description
480 !% Like <tt>XYZVelocities</tt> but in PDB format, as in <tt>PDBCoordinates</tt>.
481 !%End
482
483 !%Variable Velocities
484 !%Type block
485 !%Section System::Velocities
486 !%Description
487 !% If <tt>XYZVelocities</tt>, <tt>PDBVelocities</tt>, and <tt>XSFVelocities</tt>
488 !% are not present, <tt>Octopus</tt> will try to fetch the initial
489 !% atomic velocities from this block. If this block is not present, <tt>Octopus</tt>
490 !% will set the initial velocities to zero. The format of this block can be
491 !% illustrated by this example:
492 !%
493 !% <tt>%Velocities
494 !% <br>&nbsp;&nbsp;'C' | -1.7 | 0.0 | 0.0
495 !% <br>&nbsp;&nbsp;'O' | &nbsp;1.7 | 0.0 | 0.0
496 !% <br>%</tt>
497 !%
498 !% It describes one carbon and one oxygen moving at the relative
499 !% velocity of 3.4 velocity units.
500 !%
501 !% Note: It is important for the velocities to maintain the ordering
502 !% in which the atoms were defined in the coordinates specifications.
503 !%End
504
505 call read_coords_init(xyz)
506 call read_coords_read('Velocities', xyz, ions%space, namespace)
507 if (xyz%source /= read_coords_err) then
508
509 have_velocities = .true.
510
511 if (ions%natoms /= xyz%n) then
512 write(message(1), '(a,i4,a,i4)') 'I need exactly ', ions%natoms, ' velocities, but I found ', xyz%n
513 call messages_fatal(1, namespace=namespace)
514 end if
515
516 ! copy information and adjust units
517 do i = 1, ions%natoms
518 ions%vel(:, i) = units_to_atomic(units_inp%velocity/units_inp%length, xyz%atom(i)%x(1:ions%space%dim))
519 end do
520
521 call read_coords_end(xyz)
522
523 else
524 ions%vel = m_zero
525 end if
526 end if
527
528 call ions%update_kinetic_energy()
529
530 !%Variable MoveIons
531 !%Type logical
532 !%Section Time-Dependent::Propagation
533 !%Description
534 !% This variable controls whether atoms are moved during a time
535 !% propagation run. The default is yes when the ion velocity is
536 !% set explicitly or implicitly, otherwise is no.
537 !%End
538 call parse_variable(namespace, 'MoveIons', have_velocities, this%move_ions)
539 call messages_print_var_value('MoveIons', this%move_ions, namespace=namespace)
540
541 if (this%move_ions .and. ions%space%periodic_dim == 1) then
542 call messages_input_error(namespace, 'MoveIons', &
543 'Moving ions for a 1D periodic system is not allowed, as forces are incorrect.')
544 end if
545
546 if (this%ions_move()) then
547 safe_allocate(this%oldforce(1:ions%space%dim, 1:ions%natoms))
548 end if
549
550 if (ions%space%is_periodic()) then
551 !%Variable CellDynamics
552 !%Type logical
553 !%Section Time-Dependent::Propagation
554 !%Description
555 !% This variable controls whether the cell relaxation is done during a time
556 !% propagation run. The default is no.
557 !% This is done based on the Parrinello-Rahman equation of motion of the cell,
558 !% see Parrinello and Rahman, J. Appl. Pys. 52, 7182 (1981).
559 !% Only for periodic systems.
560 !%End
561 call parse_variable(namespace, 'CellDynamics', .false., this%relax_cell)
562 call messages_print_var_value('CellDynamics', this%relax_cell, namespace=namespace)
563
564 if (this%cell_relax()) then
565 if (accel_is_enabled()) then
566 message(1) = "Cell dynamics not supported on GPUs."
567 call messages_fatal(1, namespace=namespace)
568 end if
569
570 periodic_dim = ions%space%periodic_dim
571 ncomp = periodic_dim * periodic_dim
572 safe_allocate(this%cell_force(1:ncomp))
573 this%cell_force = m_zero
574 safe_allocate(this%old_cell_force(1:ncomp))
575 this%old_cell_force = m_zero
576 safe_allocate(this%cell_vel(1:ncomp))
577 this%cell_vel = m_zero
578 ! We start from identity for the strain
579 safe_allocate(this%strain(1:ncomp))
580 this%strain = m_zero
581 ncomp = 1
582 do i = 1, periodic_dim
583 do j = i, periodic_dim
584 if(i == j) this%strain(ncomp) = m_one
585 ncomp = ncomp + 1
586 end do
587 end do
588
589 ! As we work with "reduced" quantities, the initial lattice vectors are used as a reference
590 ! and the strain is propagated from this initial reference
591 safe_allocate(this%initial_rlattice(1:periodic_dim, 1:periodic_dim))
592 this%initial_rlattice(1:periodic_dim, 1:periodic_dim) = ions%latt%rlattice(1:periodic_dim, 1:periodic_dim)
593 end if
594
595 !%Variable HydrostaticPressure
596 !%Type float
597 !%Default 0.0
598 !%Section Time-Dependent::Propagation
599 !%Description
600 !% Geometry optimization and molecular dynamics can be performed in presence of an external
601 !% hydrostatic pressure. This variable allows to set this value.
602 !% Only for periodic systems.
603 !%End
604 call parse_variable(namespace, 'HydrostaticPressure', m_zero, this%pressure)
605 else
606 this%relax_cell = .false.
607 end if
608
609
610 pop_sub(ion_dynamics_init)
611 end subroutine ion_dynamics_init
612
613
614 ! ---------------------------------------------------------
615 subroutine ion_dynamics_end(this)
616 type(ion_dynamics_t), intent(inout) :: this
617
618 push_sub(ion_dynamics_end)
619 safe_deallocate_a(this%oldforce)
620
621 if (this%thermostat /= thermo_none) then
622 call tdf_end(this%temperature_function)
623 end if
624
625 if (this%drive_ions .and. allocated(this%td_displacements)) then
626 if (any(this%td_displacements(1:this%ions_t0%natoms)%move)) then
627 ! ions end cannot be called here, otherwise the species are destroyed twice
628 ! call ions_end(this%ions_t0)
629 !AFE_DEALLOCATE_P(this%ions_t0)
630 end if
631 safe_deallocate_a(this%td_displacements)
632 end if
633
634 safe_deallocate_a(this%cell_force)
635 safe_deallocate_a(this%old_cell_force)
636 safe_deallocate_a(this%cell_vel)
637 safe_deallocate_a(this%initial_rlattice)
638
639 pop_sub(ion_dynamics_end)
640 end subroutine ion_dynamics_end
641
642
643 ! ---------------------------------------------------------
645 subroutine ion_dynamics_propagate(this, ions, time, dt, namespace)
646 type(ion_dynamics_t), intent(inout) :: this
647 type(ions_t), intent(inout) :: ions
648 real(real64), intent(in) :: time
649 real(real64), intent(in) :: dt
650 type(namespace_t), intent(in) :: namespace
651
652 integer :: iatom
653
654 push_sub(ion_dynamics_propagate)
655
656 this%dt = dt
657
658 if (this%drive_ions) then
659
660 call ion_dynamics_propagate_driven_ions(this, ions, time, dt)
661
662 else
663
664 ! Update the thermostat temperature
665 call ion_dynamics_update_temperature(this, time, namespace)
666
667 ! Conversion to reduced coordinates
668 do iatom = 1, ions%natoms
669 ions%pos(:, iatom) = ions%latt%cart_to_red(ions%pos(:, iatom))
670 ions%vel(:, iatom) = ions%latt%cart_to_red(ions%vel(:, iatom))
671 ions%tot_force(:, iatom) = ions%latt%cart_to_red(ions%tot_force(:, iatom))
672 end do
673
674 ! Get the new reduced coordinates
675 if (this%ions_move()) then
676 call ion_dynamics_propagate_ions(this, ions, time, dt)
677 end if
678
679 ! Updating the lattice vectors
680 if (this%cell_relax()) then
681 call ion_dynamics_propagate_cell(this, ions, time, dt, namespace)
682 end if
683
684 ! Get the new Cartesian coordinates
685 do iatom = 1, ions%natoms
686 ions%pos(:, iatom) = ions%latt%red_to_cart(ions%pos(:, iatom))
687 ions%vel(:, iatom) = ions%latt%red_to_cart(ions%vel(:, iatom))
688 if (allocated(this%oldforce)) then
689 this%oldforce(:, iatom) = ions%latt%red_to_cart(this%oldforce(:, iatom))
690 end if
691 ions%tot_force(:, iatom) = ions%latt%red_to_cart(ions%tot_force(:, iatom))
692 end do
693
694 end if
695
696 call ions%fold_atoms_into_cell()
697
699 end subroutine ion_dynamics_propagate
700
701
702 ! ---------------------------------------------------------
704 subroutine ion_dynamics_update_temperature(this, time, namespace)
705 type(ion_dynamics_t), intent(inout) :: this
706 real(real64), intent(in) :: time
707 type(namespace_t), intent(in) :: namespace
710
711 ! get the temperature from the tdfunction for the current time
712 if (this%thermostat /= thermo_none) then
713 this%current_temperature = units_to_atomic(unit_kelvin, tdf(this%temperature_function, time))
714
715 if (this%current_temperature < m_zero) then
716 write(message(1), '(a, f10.3, 3a, f10.3, 3a)') &
717 "Negative temperature (", &
718 units_from_atomic(unit_kelvin, this%current_temperature), " ", units_abbrev(unit_kelvin), &
719 ") at time ", &
720 units_from_atomic(units_out%time, time), " ", trim(units_abbrev(units_out%time)), "."
721 call messages_fatal(1, namespace=namespace)
722 end if
723 else
724 this%current_temperature = m_zero
725 end if
726
729
733 subroutine ion_dynamics_propagate_driven_ions(this, ions, time, dt)
734 type(ion_dynamics_t), intent(inout) :: this
735 type(ions_t), intent(inout) :: ions
736 real(real64), intent(in) :: time
737 real(real64), intent(in) :: dt
739 integer :: iatom
740 real(real64) :: dr(3)
741
743
744 assert(this%drive_ions)
745
746 do iatom = 1, ions%natoms
747 if (ions%fixed(iatom)) cycle
748
749 if (this%constant_velocity) then
750 ions%pos(:, iatom) = ions%pos(:, iatom) + dt*ions%vel(:, iatom)
751
752 else if (allocated(this%td_displacements)) then
753
754 if (this%td_displacements(iatom)%move) then
755 dr(1:3)=(/ real(tdf(this%td_displacements(iatom)%fx, time), real64), &
756 real(tdf(this%td_displacements(iatom)%fy, time), real64), &
757 real(tdf(this%td_displacements(iatom)%fz, time), real64) /)
758
759 ions%pos(:, iatom) = this%ions_t0%pos(:, iatom) + dr(1:ions%space%dim)
760 end if
761
762 end if
763 end do
764
767
768 ! ---------------------------------------------------------
773 subroutine ion_dynamics_propagate_ions(this, ions, time, dt)
774 type(ion_dynamics_t), intent(inout) :: this
775 type(ions_t), intent(inout) :: ions
776 real(real64), intent(in) :: time
777 real(real64), intent(in) :: dt
778
779 integer :: iatom
780
782
783 assert(.not. this%drive_ions)
784
785 if (this%thermostat /= thermo_nh) then
786 ! integrate using verlet
787 do iatom = 1, ions%natoms
788 if (ions%fixed(iatom)) cycle
789
790 ions%pos(:, iatom) = ions%pos(:, iatom) + dt*ions%vel(:, iatom) + &
791 m_half*dt**2 / ions%mass(iatom) * ions%tot_force(:, iatom)
792
793 this%oldforce(:, iatom) = ions%tot_force(:, iatom)
794 end do
795
796 else
797 ! for the Nose-Hoover thermostat we use a special integrator
798
799 ! The implementation of the Nose-Hoover thermostat is based on
800 ! Understanding Molecular Simulations by Frenkel and Smit,
801 ! Appendix E, page 540-542.
802
803 call nh_chain(this, ions)
804
805 do iatom = 1, ions%natoms
806 if (ions%fixed(iatom)) cycle
807
808 ions%pos(:, iatom) = ions%pos(:, iatom) + m_half*dt*ions%vel(:, iatom)
809 end do
810
811 end if
812
814 end subroutine ion_dynamics_propagate_ions
815
816 ! ---------------------------------------------------------
818 subroutine ion_dynamics_propagate_cell(this, ions, time, dt, namespace)
819 type(ion_dynamics_t), intent(inout) :: this
820 type(ions_t), intent(inout) :: ions
821 real(real64), intent(in) :: time
822 real(real64), intent(in) :: dt
823 type(namespace_t), intent(in) :: namespace
824
825 integer :: idir, jdir, comp
826 real(real64) :: rlattice_change(ions%space%periodic_dim*ions%space%periodic_dim)
827
829
830 rlattice_change = dt * this%cell_vel + m_half*dt**2 * this%cell_force
831
832 write(message(1),'(a,3a,a)') ' Cell force [', trim(units_abbrev(units_out%energy/units_out%length**ions%space%dim)), ']'
833 do idir = 1, ions%space%periodic_dim
834 write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%energy/units_out%length**ions%space%dim, &
835 this%cell_force(jdir + (idir-1)*ions%space%periodic_dim)), &
836 jdir = 1, ions%space%periodic_dim)
837 end do
838 call messages_info(1+ions%space%periodic_dim, namespace=ions%namespace)
839
840 write(message(1),'(a,3a,a)') ' Cell vel [', &
841 trim(units_abbrev(units_out%energy/units_out%length**ions%space%dim*units_out%time)), ']'
842 do idir = 1, ions%space%periodic_dim
843 write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%energy/units_out%length**ions%space%dim*units_out%time, &
844 this%cell_vel(jdir+ (idir-1)*ions%space%periodic_dim)), &
845 jdir = 1, ions%space%periodic_dim)
846 end do
847 call messages_info(1+ions%space%periodic_dim, namespace=ions%namespace)
848
849
850 comp = 1
851 do idir = 1, ions%space%periodic_dim
852 do jdir = 1, ions%space%periodic_dim
853 ions%latt%rlattice(idir, jdir) = ions%latt%rlattice(idir, jdir) + rlattice_change(comp)
854 comp = comp + 1
855 end do
856 end do
857
858 this%old_cell_force = this%cell_force
859
860 if (associated(this%symm)) then
861 call this%symm%symmetrize_lattice_vectors(ions%space%periodic_dim, &
862 this%initial_rlattice, ions%latt%rlattice, this%symmetrize)
863 end if
864 call ions%update_lattice_vectors(ions%latt, this%symmetrize)
865
867 end subroutine ion_dynamics_propagate_cell
868
869
870 ! ---------------------------------------------------------
871 subroutine nh_chain(this, ions)
872 type(ion_dynamics_t), intent(inout) :: this
873 type(ions_t), intent(inout) :: ions
874
875 real(real64) :: g1, g2, ss, uk, dt, temp
876
877 push_sub(nh_chain)
878
879 dt = this%dt
880
881 call ions%update_kinetic_energy()
882 uk = ions%kinetic_energy
883
884 temp = this%current_temperature
885
886 g2 = (this%nh(1)%mass*this%nh(1)%vel**2 - temp)/this%nh(2)%mass
887 this%nh(2)%vel = this%nh(2)%vel + g2*dt/m_four
888 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
889
890 g1 = (m_two*uk - m_three*ions%natoms*temp)/this%nh(1)%mass
891 this%nh(1)%vel = this%nh(1)%vel + g1*dt/m_four
892 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
893 this%nh(1)%pos = this%nh(1)%pos + this%nh(1)%vel*dt/m_two
894 this%nh(2)%pos = this%nh(2)%pos + this%nh(2)%vel*dt/m_two
895
896 ss = exp(-this%nh(1)%vel*dt/m_two)
897
898 ions%vel = ss*ions%vel
899
900 uk = uk*ss**2
901
902 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
903 g1 = (m_two*uk - m_three*ions%natoms*temp)/this%nh(1)%mass
904 this%nh(1)%vel = this%nh(1)%vel + g1*dt/m_four
905 this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
906
907 g2 = (this%nh(1)%mass*this%nh(1)%vel**2 - temp)/this%nh(2)%mass
908 this%nh(2)%vel = this%nh(2)%vel + g2*dt/m_four
909
910 pop_sub(nh_chain)
911 end subroutine nh_chain
912
913
914 ! ---------------------------------------------------------
915 subroutine ion_dynamics_propagate_vel(this, ions, atoms_moved)
916 type(ion_dynamics_t), intent(inout) :: this
917 type(ions_t), intent(inout) :: ions
918 logical, optional, intent(out) :: atoms_moved
919
920 integer :: iatom
921 real(real64) :: scal
922
923 if (.not. ion_dynamics_ions_move(this)) return
924 if (this%drive_ions) return
925
927
928 if (present(atoms_moved)) atoms_moved = this%thermostat == thermo_nh
929
930 if (this%thermostat /= thermo_nh) then
931 ! velocity verlet
932
933 do iatom = 1, ions%natoms
934 if (ions%fixed(iatom)) cycle
935
936 ions%vel(:, iatom) = ions%vel(:, iatom) &
937 + this%dt/ions%mass(iatom) * m_half * (this%oldforce(:, iatom) + &
938 ions%tot_force(:, iatom))
939
940 end do
941
942 else
943 ! the nose-hoover integration
944 do iatom = 1, ions%natoms
945 ions%vel(:, iatom) = ions%vel(:, iatom) + this%dt*ions%tot_force(:, iatom) / ions%mass(iatom)
946 ions%pos(:, iatom) = ions%pos(:, iatom) + m_half*this%dt*ions%vel(:, iatom)
947 end do
948
949 call nh_chain(this, ions)
950
951 end if
952
953 if (this%thermostat == thermo_scal) then
954 scal = sqrt(this%current_temperature/ion_dynamics_temperature(ions))
955
956 ions%vel = scal*ions%vel
957 end if
958
959 if (this%cell_relax()) then
960 this%cell_vel = this%cell_vel + this%dt * m_half * (this%old_cell_force + this%cell_force)
961 end if
962
965
966
967 ! ---------------------------------------------------------
969 subroutine ion_dynamics_verlet_step1(ions, q, v, fold, dt)
970 type(ions_t), intent(in) :: ions
971 real(real64), intent(inout) :: q(:, :)
972 real(real64), intent(inout) :: v(:, :)
973 real(real64), intent(in) :: fold(:, :)
974 real(real64), intent(in) :: dt
975
976 integer :: iatom
977
979
980 ! First transform momenta to velocities
981 do iatom = 1, ions%natoms
982 v(iatom, 1:ions%space%dim) = v(iatom, 1:ions%space%dim) / ions%mass(iatom)
983 end do
984
985 ! integrate using verlet
986 do iatom = 1, ions%natoms
987 if (ions%fixed(iatom)) cycle
988 q(iatom, 1:ions%space%dim) = q(iatom, 1:ions%space%dim) + dt * v(iatom, 1:ions%space%dim) + &
989 m_half*dt**2 / ions%mass(iatom) * fold(iatom, 1:ions%space%dim)
990 end do
991
992 ! And back to momenta.
993 do iatom = 1, ions%natoms
994 v(iatom, 1:ions%space%dim) = ions%mass(iatom) * v(iatom, 1:ions%space%dim)
995 end do
996
998 end subroutine ion_dynamics_verlet_step1
999
1000
1001
1002 ! ---------------------------------------------------------
1004 subroutine ion_dynamics_verlet_step2(ions, v, fold, fnew, dt)
1005 type(ions_t), intent(in) :: ions
1006 real(real64), intent(inout) :: v(:, :)
1007 real(real64), intent(in) :: fold(:, :)
1008 real(real64), intent(in) :: fnew(:, :)
1009 real(real64), intent(in) :: dt
1010
1011 integer :: iatom
1012
1014
1015 ! First transform momenta to velocities
1016 do iatom = 1, ions%natoms
1017 v(iatom, 1:ions%space%dim) = v(iatom, 1:ions%space%dim) / ions%mass(iatom)
1018 end do
1019
1020 ! velocity verlet
1021 do iatom = 1, ions%natoms
1022 if (ions%fixed(iatom)) cycle
1023 v(iatom, 1:ions%space%dim) = v(iatom, 1:ions%space%dim) &
1024 + dt / ions%mass(iatom) * m_half * (fold(iatom, 1:ions%space%dim) + fnew(iatom, 1:ions%space%dim))
1025 end do
1026
1027 ! And back to momenta.
1028 do iatom = 1, ions%natoms
1029 v(iatom, 1:ions%space%dim) = ions%mass(iatom) * v(iatom, 1:ions%space%dim)
1030 end do
1031
1033 end subroutine ion_dynamics_verlet_step2
1034
1035
1036 ! ---------------------------------------------------------
1037 subroutine ion_dynamics_save_state(this, ions, state)
1038 type(ion_dynamics_t), intent(in) :: this
1039 type(ions_t), intent(in) :: ions
1040 type(ion_state_t), intent(out) :: state
1041
1042 if (.not. this%ions_move()) return
1043
1044 push_sub(ion_dynamics_save_state)
1045
1046 safe_allocate(state%pos(1:ions%space%dim, 1:ions%natoms))
1047 safe_allocate(state%vel(1:ions%space%dim, 1:ions%natoms))
1048
1049 state%pos = ions%pos
1050 state%vel = ions%vel
1051
1052 if (this%thermostat == thermo_nh) then
1053 safe_allocate(state%old_pos(1:ions%space%dim, 1:ions%natoms))
1054 state%old_pos(1:ions%space%dim, 1:ions%natoms) = this%old_pos(1:ions%space%dim, 1:ions%natoms)
1055 state%nh(1:2)%pos = this%nh(1:2)%pos
1056 state%nh(1:2)%vel = this%nh(1:2)%vel
1057 end if
1058
1060 end subroutine ion_dynamics_save_state
1061
1063 ! ---------------------------------------------------------
1064 subroutine ion_dynamics_restore_state(this, ions, state)
1065 type(ion_dynamics_t), intent(inout) :: this
1066 type(ions_t), intent(inout) :: ions
1067 type(ion_state_t), intent(inout) :: state
1068
1069 assert(.not. this%cell_relax())
1070 if (.not. this%ions_move()) return
1071
1073
1074 ions%pos = state%pos
1075 ions%vel = state%vel
1076
1077 if (this%thermostat == thermo_nh) then
1078 this%old_pos(1:ions%space%dim, 1:ions%natoms) = state%old_pos(1:ions%space%dim, 1:ions%natoms)
1079 this%nh(1:2)%pos = state%nh(1:2)%pos
1080 this%nh(1:2)%vel = state%nh(1:2)%vel
1081 safe_deallocate_a(state%old_pos)
1082 end if
1083
1084 safe_deallocate_a(state%pos)
1085 safe_deallocate_a(state%vel)
1086
1088 end subroutine ion_dynamics_restore_state
1089
1090
1091 ! ---------------------------------------------------------
1092 logical pure function ion_dynamics_ions_move(this) result(ions_move)
1093 class(ion_dynamics_t), intent(in) :: this
1094
1095 ions_move = this%move_ions
1096
1098
1099
1100 ! ---------------------------------------------------------
1102 logical pure function ion_dynamics_drive_ions(this) result(drive_ions)
1103 type(ion_dynamics_t), intent(in) :: this
1104
1105 drive_ions = this%drive_ions
1106
1107 end function ion_dynamics_drive_ions
1108
1109 ! ---------------------------------------------------------
1111 logical pure function ion_dynamics_cell_relax(this) result(cell_dynamics)
1112 class(ion_dynamics_t), intent(in) :: this
1113
1114 cell_dynamics = this%relax_cell
1115
1116 end function ion_dynamics_cell_relax
1117
1118 ! ---------------------------------------------------------
1120 logical pure function ion_dynamics_is_active(this) result(is_active)
1121 class(ion_dynamics_t), intent(in) :: this
1122
1123 is_active = this%relax_cell .or. this%move_ions
1124
1125 end function ion_dynamics_is_active
1126
1127
1128 ! ---------------------------------------------------------
1130 real(real64) function ion_dynamics_temperature(ions) result(temperature)
1131 type(ions_t), intent(in) :: ions
1132
1133 integer :: iatom
1134 real(real64) :: kinetic_energy
1135
1136 kinetic_energy = m_zero
1137 do iatom = 1, ions%natoms
1138 kinetic_energy = kinetic_energy + &
1139 m_half * ions%mass(iatom) * sum(ions%vel(:, iatom)**2)
1140 end do
1141 temperature = m_two/m_three*kinetic_energy/ions%natoms
1142
1143 end function ion_dynamics_temperature
1144
1145
1146 ! ---------------------------------------------------------
1148 logical function ion_dynamics_freeze(this) result(freeze)
1149 type(ion_dynamics_t), intent(inout) :: this
1150 if (this%move_ions) then
1151 this%move_ions = .false.
1152 freeze = .true.
1153 else
1154 freeze = .false.
1155 end if
1156 end function ion_dynamics_freeze
1158
1159 ! ---------------------------------------------------------
1161 subroutine ion_dynamics_unfreeze(this)
1162 type(ion_dynamics_t), intent(inout) :: this
1163 this%move_ions = .true.
1164 end subroutine ion_dynamics_unfreeze
1165
1166 ! ---------------------------------------------------------
1167 subroutine ion_dynamics_dump(this, restart, ierr)
1168 type(ion_dynamics_t), intent(in) :: this
1169 type(restart_t), intent(in) :: restart
1170 integer, intent(out) :: ierr
1171
1172 push_sub(ion_dynamics_dump)
1173
1174 if (allocated(this%oldforce)) then
1175 call drestart_write_binary(restart, "ion_dynamics_oldforce", size(this%oldforce), &
1176 this%oldforce, ierr)
1177 end if
1178
1179 if( allocated(this%old_cell_force)) then
1180 call drestart_write_binary(restart, "ion_dynamics_old_cell_force", size(this%old_cell_force), &
1181 this%old_cell_force, ierr)
1182 end if
1183
1184 pop_sub(ion_dynamics_dump)
1185 end subroutine ion_dynamics_dump
1186
1187 ! ---------------------------------------------------------
1188 subroutine ion_dynamics_load(this, restart, ierr)
1189 type(ion_dynamics_t), intent(inout) :: this
1190 type(restart_t), intent(in) :: restart
1191 integer, intent(out) :: ierr
1192
1193 push_sub(ion_dynamics_load)
1194
1195 if (allocated(this%oldforce)) then
1196 call drestart_read_binary(restart, "ion_dynamics_oldforce", size(this%oldforce), &
1197 this%oldforce, ierr)
1198 end if
1199
1200 if( allocated(this%old_cell_force)) then
1201 call drestart_read_binary(restart, "ion_dynamics_old_cell_force", size(this%old_cell_force), &
1202 this%old_cell_force, ierr)
1203 end if
1205 pop_sub(ion_dynamics_load)
1206 end subroutine ion_dynamics_load
1207
1208 !----------------------------------------------------------
1210 subroutine ion_dynamics_update_stress(this, space, stress, rlattice, rcell_volume)
1211 class(ion_dynamics_t), intent(inout) :: this
1212 class(space_t), intent(in) :: space
1213 real(real64), intent(in) :: stress(3,3)
1214 real(real64), intent(in) :: rlattice(:,:)
1215 real(real64), intent(in) :: rcell_volume
1216
1217 integer :: idir, jdir, comp
1218 real(real64) :: inv_latt(space%periodic_dim, space%periodic_dim), tmp_stress(space%periodic_dim, space%periodic_dim)
1219 real(real64) :: cell_force(space%periodic_dim, space%periodic_dim)
1220
1222
1223 ! Get the inverse lattice
1224 inv_latt = rlattice(1:space%periodic_dim, 1:space%periodic_dim)
1225 call lalg_inverse(space%periodic_dim, inv_latt, 'dir')
1226
1227 tmp_stress = -stress(1:space%periodic_dim, 1:space%periodic_dim)
1228 do idir = 1, space%periodic_dim
1229 tmp_stress(idir, idir) = tmp_stress(idir, idir) - this%pressure/space%periodic_dim
1230 end do
1231 cell_force = matmul(tmp_stress, transpose(inv_latt)) * rcell_volume
1232
1233 comp = 1
1234 do idir = 1, space%periodic_dim
1235 do jdir = 1, space%periodic_dim
1236 this%cell_force(comp) = cell_force(idir, jdir)
1237 comp = comp + 1
1238 end do
1239 end do
1240
1241 if (debug%info) then
1242 write(message(1),'(a,3a,a)') ' Stress tensor [', trim(units_abbrev(units_out%energy/units_out%length**space%dim)), ']'
1243 do idir = 1, space%periodic_dim
1244 write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%energy/units_out%length**space%dim, stress(jdir, idir)), &
1245 jdir = 1, space%periodic_dim)
1246 end do
1247 call messages_info(1+space%periodic_dim, namespace=global_namespace)
1248 end if
1249
1251 end subroutine ion_dynamics_update_stress
1252
1253 !----------------------------------------------------------
1254 subroutine electrons_lattice_vectors_update(namespace, gr, space, psolver, kpoints, mc, qtot, new_latt)
1255 type(namespace_t), intent(in) :: namespace
1256 type(grid_t), intent(inout) :: gr
1257 class(space_t), intent(in) :: space
1258 type(poisson_t), intent(inout) :: psolver
1259 type(kpoints_t), intent(inout) :: kpoints
1260 type(multicomm_t), intent(in) :: mc
1261 real(real64), intent(in) :: qtot
1262 type(lattice_vectors_t), intent(in) :: new_latt
1263
1264 integer :: idir
1265 real(real64) :: length(1:space%dim)
1266
1268
1269 ! Regenerate the box
1270 select type(box => gr%box)
1271 type is (box_parallelepiped_t)
1272 do idir = 1, space%dim
1273 length(idir) = norm2(new_latt%rlattice(1:space%dim, idir))
1274 end do
1275 call box%regenerate(space%dim, new_latt%rlattice, length, namespace)
1276 class default
1277 call messages_not_implemented("Grid regeneration for non-parallelepiped boxes", namespace=namespace)
1278 end select
1279
1280 call grid_lattice_vectors_update(gr, space, namespace, mc, new_latt)
1282 !Initialize Poisson solvers
1283 call poisson_end(psolver)
1284 call poisson_init(psolver, namespace, space, gr%der, mc, gr%stencil, qtot, verbose=.false.)
1285
1286 call kpoints_lattice_vectors_update(kpoints, new_latt)
1287
1289
1291
1292
1293end module ion_dynamics_oct_m
1294
1295!! Local Variables:
1296!! mode: f90
1297!! coding: utf-8
1298!! End:
Functions to generate random numbers.
Definition: loct_math.F90:329
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:180
double exp(double __x) __attribute__((__nothrow__
pure logical function, public accel_is_enabled()
Definition: accel.F90:402
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public ion_dynamics_verlet_step2(ions, v, fold, fnew, dt)
A bare verlet integrator.
subroutine ion_dynamics_update_temperature(this, time, namespace)
Update the temperature of the ions in case of a thermostat.
subroutine, public ion_dynamics_dump(this, restart, ierr)
subroutine, public ion_dynamics_verlet_step1(ions, q, v, fold, dt)
A bare verlet integrator.
logical pure function ion_dynamics_cell_relax(this)
Is the cell dynamics activated or not.
subroutine nh_chain(this, ions)
subroutine, public ion_dynamics_restore_state(this, ions, state)
subroutine, public ion_dynamics_propagate(this, ions, time, dt, namespace)
Interface for the ion/cell dynamics.
subroutine ion_dynamics_propagate_cell(this, ions, time, dt, namespace)
Time-evolution of the lattice vectors.
subroutine, public ion_dynamics_save_state(this, ions, state)
subroutine, public ion_dynamics_unfreeze(this)
Unfreezes the ionic movement.
subroutine ion_dynamics_propagate_ions(this, ions, time, dt)
Time evolution of the ions.
subroutine ion_dynamics_update_stress(this, space, stress, rlattice, rcell_volume)
Updates the stress tensor for the ion dynamics.
subroutine, public ion_dynamics_propagate_vel(this, ions, atoms_moved)
subroutine, public ion_dynamics_load(this, restart, ierr)
subroutine, public ion_dynamics_init(this, namespace, ions, symmetrize, symm)
integer, parameter thermo_scal
subroutine, public electrons_lattice_vectors_update(namespace, gr, space, psolver, kpoints, mc, qtot, new_latt)
subroutine, public ion_dynamics_end(this)
subroutine ion_dynamics_propagate_driven_ions(this, ions, time, dt)
Move ions following a driven motion.
integer, parameter thermo_nh
logical pure function ion_dynamics_ions_move(this)
logical function, public ion_dynamics_freeze(this)
Freezes the ionic movement.
logical pure function ion_dynamics_is_active(this)
Is the cell dynamics activated or not.
logical pure function, public ion_dynamics_drive_ions(this)
Is the ion dynamics activated or not.
real(real64) function, public ion_dynamics_temperature(ions)
This function returns the ionic temperature in energy units.
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:434
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:270
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:502
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
integer, parameter, public read_coords_err
for read_coords_info::file_type
subroutine, public read_coords_init(gf)
subroutine, public read_coords_end(gf)
subroutine, public read_coords_read(what, gf, space, namespace)
subroutine, public tdf_end(f)
Definition: tdfunction.F90:974
subroutine, public tdf_read(f, namespace, function_name, ierr)
This function initializes "f" from the TDFunctions block.
Definition: tdfunction.F90:218
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:223
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
type(unit_t), public unit_kelvin
For converting energies into temperatures.
type(unit_system_t), public units_inp
the units systems for reading and writing
int true(void)