Octopus
dftb.F90
Go to the documentation of this file.
1!! Copyright (C) 2020 F. Bonafé, H. Appel
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 dftb_oct_m
23 use debug_oct_m
24#ifdef HAVE_DFTBPLUS
25 use dftbplus
26#endif
27 use global_oct_m
31 use io_oct_m
33 use ions_oct_m
34 use, intrinsic :: iso_fortran_env
35 use lasers_oct_m
37 use mpi_oct_m
40 use parser_oct_m
46 use system_oct_m
48 use unit_oct_m
51
52 implicit none
53
54 private
55 public :: &
56 dftb_t, &
58
59 integer, parameter :: &
60 OUTPUT_COORDINATES = 1, &
62
65 type, extends(system_t) :: dftb_t
66 integer :: n_atom
67 real(real64), allocatable :: coords(:,:), gradients(:,:)
68 real(real64), allocatable :: acc(:,:)
69 real(real64), allocatable :: tot_force(:,:)
70 real(real64), allocatable :: vel(:,:)
71 real(real64), allocatable :: prev_tot_force(:,:)
72 integer, allocatable :: species(:)
73 integer :: dynamics
74 real(real64), allocatable :: mass(:)
75 real(real64), allocatable :: atom_charges(:,:)
76 character(len=LABEL_LEN), allocatable :: labels(:)
77 real(real64), allocatable :: prev_acc(:,:,:)
78 real(real64) :: scc_tolerance
79 class(ions_t), pointer :: ions => null()
80 type(c_ptr) :: output_handle(2)
81 type(ion_dynamics_t) :: ions_dyn
82 class(lasers_t), pointer :: lasers => null()
83 logical :: laser_field
84 real(real64) :: field(3)
85 real(real64) :: energy
86#ifdef HAVE_DFTBPLUS
87 type(TDftbPlus) :: dftbp
88#endif
89 contains
90 procedure :: init_interaction => dftb_init_interaction
91 procedure :: initialize => dftb_initialize
92 procedure :: do_algorithmic_operation => dftb_do_algorithmic_operation
93 procedure :: output_start => dftb_output_start
94 procedure :: output_write => dftb_output_write
95 procedure :: output_finish => dftb_output_finish
96 procedure :: is_tolerance_reached => dftb_is_tolerance_reached
97 procedure :: copy_quantities_to_interaction => dftb_copy_quantities_to_interaction
98 procedure :: init_interaction_as_partner => dftb_init_interaction_as_partner
99 procedure :: update_interactions_start => dftb_update_interactions_start
100 procedure :: restart_write_data => dftb_restart_write_data
101 procedure :: restart_read_data => dftb_restart_read_data
102 procedure :: update_kinetic_energy => dftb_update_kinetic_energy
103 final :: dftb_finalize
104 end type dftb_t
105
106 interface dftb_t
107 procedure dftb_constructor
108 end interface dftb_t
109
111 integer, parameter :: &
112 EHRENFEST = 1, &
113 bo = 2
115contains
116
117 ! ---------------------------------------------------------
122 function dftb_constructor(namespace) result(sys)
123 class(dftb_t), pointer :: sys
124 type(namespace_t), intent(in) :: namespace
125
126 push_sub(dftb_constructor)
127
128#ifndef HAVE_DFTBPLUS
129 message(1) = "DFTB+ system not available. This feature requires compiling Octopus with the DFTB+ library"
130 call messages_fatal(1, namespace=namespace)
131#endif
132
133 allocate(sys)
134
135 call dftb_init(sys, namespace)
136
137 pop_sub(dftb_constructor)
138 end function dftb_constructor
139
140 ! ---------------------------------------------------------
145 ! ---------------------------------------------------------
146 subroutine dftb_init(this, namespace)
147 class(dftb_t), target, intent(inout) :: this
148 type(namespace_t), intent(in) :: namespace
149
150 integer :: ii, jj, ispec
151 character(len=MAX_PATH_LEN) :: slako_dir
152 character(len=1), allocatable :: max_ang_mom(:)
153 character(len=LABEL_LEN) :: this_max_ang_mom, this_label
154 integer :: n_maxang_block
155 type(block_t) :: blk
156 real(real64) :: initial_temp
157
158#ifdef HAVE_DFTBPLUS
159 type(tdftbplusinput) :: input
160 type(fnode), pointer :: proot, pgeo, pham, pdftb, pmaxang, pslakos, ptype2files, panalysis
161 type(fnode), pointer :: pparseropts
162 type(fnode), pointer :: pelecdyn, pperturb, plaser
163#endif
165 push_sub(dftb_init)
167 this%namespace = namespace
169 ! Currently this system does not support any interaction
170 allocate(this%supported_interactions(0))
171 allocate(this%supported_interactions_as_partner(0))
173 ! Quantities updated by the algorithm
174 call this%quantities%add(quantity_t("position", updated_on_demand = .false.))
175 call this%quantities%add(quantity_t("velocity", updated_on_demand = .false.))
177 call messages_print_with_emphasis(msg="DFTB+ System", namespace=namespace)
179 this%ions => ions_t(namespace)
180 this%n_atom = this%ions%natoms
181 safe_allocate(this%coords(1:3, 1:this%n_atom))
182 safe_allocate(this%acc(1:3, 1:this%n_atom))
183 safe_allocate(this%vel(1:3, 1:this%n_atom))
184 safe_allocate(this%tot_force(1:3, 1:this%n_atom))
185 safe_allocate(this%prev_tot_force(1:3, 1:this%n_atom))
186 safe_allocate(this%gradients(1:3, 1:this%n_atom))
187 safe_allocate(this%species(1:this%n_atom))
188 safe_allocate(this%mass(1:this%n_atom))
189 safe_allocate(this%atom_charges(1:this%n_atom, 1))
190 safe_allocate(this%labels(1:this%ions%nspecies))
191 safe_allocate(max_ang_mom(1:this%ions%nspecies))
193 ispec = 1
194 this%species(1) = 1
195 this%labels(1) = trim(this%ions%atom(1)%label)
197 this%coords = this%ions%pos
198 ! mass is read from the default pseudopotential files
199 this%mass = this%ions%mass
200 do ii = 1, this%n_atom
201 if ((ii > 1) .and. .not. (any(this%labels(1:ispec) == this%ions%atom(ii)%label))) then
202 ispec = ispec + 1
203 this%labels(ispec) = trim(this%ions%atom(ii)%label)
204 end if
205 do jj = 1, ispec
206 if (trim(this%ions%atom(ii)%label) == trim(this%labels(jj))) then
207 this%species(ii) = jj
208 end if
209 end do
210 end do
211 this%vel = m_zero
212 this%tot_force = m_zero
213
214 !%Variable MaxAngularMomentum
215 !%Type block
216 !%Section DFTBPlusInterface
217 !%Description
218 !% Specifies the highest angular momentum for each atom type. All orbitals up
219 !% to that angular momentum will be included in the calculation.
220 !% Possible values for the angular momenta are s, p, d, f.
221 !% These are examples:
222 !%
223 !% <tt>%MaxAngularMomentum
224 !% <br>&nbsp;&nbsp;'O' | 'p'
225 !% <br>&nbsp;&nbsp;'H' | 's'
226 !% <br>%</tt>
227 !%End
228 n_maxang_block = 0
229 if (parse_block(namespace, 'MaxAngularMomentum', blk) == 0) then
230 n_maxang_block = parse_block_n(blk)
231 if (n_maxang_block /= this%ions%nspecies) then
232 call messages_input_error(namespace, "MaxAngularMomentum", "Wrong number of species.")
233 end if
234
235 do ii = 1, n_maxang_block
236 call parse_block_string(blk, ii-1, 0, this_label)
237 call parse_block_string(blk, ii-1, 1, this_max_ang_mom)
238 if (any(["s","p","d","f"] == trim(this_max_ang_mom))) then
239 call messages_input_error(namespace, "MaxAngularMomentum", "Wrong maximum angular momentum for element"//trim(this_label))
240 end if
241 do jj = 1, this%ions%nspecies
242 if (trim(adjustl(this_label)) == trim(adjustl(this%labels(jj)))) then
243 max_ang_mom(jj) = trim(adjustl(this_max_ang_mom))
244 end if
245 end do
246 end do
247 end if
248 call parse_block_end(blk)
249
250 !%Variable SlakoDir
251 !%Type string
252 !%Default "./"
253 !%Section Execution::IO
254 !%Description
255 !% Folder containing the Slako files
256 !%End
257 call parse_variable(namespace, 'SlakoDir', './', slako_dir)
258
259
260 ! Dynamics variables
261
262 call ion_dynamics_init(this%ions_dyn, namespace, this%ions)
263
264 call parse_variable(namespace, 'TDDynamics', bo, this%dynamics)
265 call messages_print_var_option('TDDynamics', this%dynamics, namespace=namespace)
266 if (this%dynamics == bo) then
267 call ion_dynamics_unfreeze(this%ions_dyn)
268 end if
269
270 !%Variable InitialIonicTemperature
271 !%Type float
272 !%Default 0.0
273 !%Section DFTBPlusInterface
274 !%Description
275 !% If this variable is present, the ions will have initial velocities
276 !% velocities to the atoms following a Boltzmann distribution with
277 !% this temperature (in Kelvin). Used only if <tt>TDDynamics = Ehrenfest</tt>
278 !% and <tt>MoveIons = yes</tt>.
279 !%End
280 call parse_variable(namespace, 'InitialIonicTemperature', m_zero, initial_temp, unit = unit_kelvin)
281
282 this%lasers => lasers_t(namespace)
283 call lasers_parse_external_fields(this%lasers)
284 if (this%lasers%no_lasers > 0) then
285 this%laser_field = .true.
286 else
287 this%laser_field = .false.
288 end if
289
290#ifdef HAVE_DFTBPLUS
291 call tdftbplus_init(this%dftbp)
292
293 call this%dftbp%getEmptyInput(input)
294 call input%getRootNode(proot)
295 call setchild(proot, "Geometry", pgeo)
296 call setchildvalue(pgeo, "Periodic", .false.)
297 call setchildvalue(pgeo, "TypeNames", this%labels(1:this%ions%nspecies))
298 call setchildvalue(pgeo, "TypesAndCoordinates", reshape(this%species, [1, size(this%species)]), this%coords)
299 call setchild(proot, "Hamiltonian", pham)
300 call setchild(pham, "Dftb", pdftb)
301 call setchildvalue(pdftb, "Scc", .true.)
302
303 !%Variable SccTolerance
304 !%Type float
305 !%Section DFTBPlusInterface
306 !%Description
307 !% Self-consistent-charges convergence tolerance. Once this
308 !% tolerance has been achieved the SCC cycle will stop.
309 !%End
310 call parse_variable(namespace, 'SccTolerance', 1e-9_real64, this%scc_tolerance)
311 call messages_print_var_value('SccTolerance', this%scc_tolerance, namespace=namespace)
312 call setchildvalue(pdftb, "SccTolerance", this%scc_tolerance)
313
314 ! sub-block inside hamiltonian for the maximum angular momenta
315 call setchild(pdftb, "MaxAngularMomentum", pmaxang)
316 ! explicitly set the maximum angular momenta for the species
317 do ii = 1, this%ions%nspecies
318 call setchildvalue(pmaxang, this%labels(ii), max_ang_mom(ii))
319 end do
320
321 ! get the SK data
322 ! You should provide the skfiles as found in the external/slakos/origin/mio-1-1/ folder. These can
323 ! be downloaded with the utils/get_opt_externals script
324 call setchild(pdftb, "SlaterKosterFiles", pslakos)
325 call setchild(pslakos, "Type2FileNames", ptype2files)
326 call setchildvalue(ptype2files, "Prefix", slako_dir)
327 call setchildvalue(ptype2files, "Separator", "-")
328 call setchildvalue(ptype2files, "Suffix", ".skf")
329
330 ! set up analysis options
331 call setchild(proot, "Analysis", panalysis)
332 call setchildvalue(panalysis, "CalculateForces", .true.)
333
334 call setchild(proot, "ParserOptions", pparseropts)
335 call setchildvalue(pparseropts, "ParserVersion", 5)
336
337 if (this%dynamics == ehrenfest) then
338 call setchild(proot, "ElectronDynamics", pelecdyn)
339 call setchildvalue(pelecdyn, "IonDynamics", ion_dynamics_ions_move(this%ions_dyn))
340 if (ion_dynamics_ions_move(this%ions_dyn)) then
341 call setchildvalue(pelecdyn, "InitialTemperature", initial_temp)
342 end if
343
344 ! initialize with wrong arguments for the moment, will be overriden later
345 call setchildvalue(pelecdyn, "Steps", 1)
346 call setchildvalue(pelecdyn, "TimeStep", m_one)
347 call setchild(pelecdyn, "Perturbation", pperturb)
348 if (this%laser_field) then
349 call setchild(pperturb, "Laser", plaser)
350 call setchildvalue(plaser, "PolarizationDirection", [ m_one , m_zero , m_zero ])
351 call setchildvalue(plaser, "LaserEnergy", m_one)
352 call setchildvalue(pelecdyn, "FieldStrength", m_one)
353 else
354 call setchild(pperturb, "None", plaser)
355 end if
356 end if
357
358 message(1) = 'Input tree in HSD format:'
359 call messages_info(1, namespace=namespace)
360 call dumphsd(input%hsdTree, stdout)
361
362 ! initialise the DFTB+ calculator
363 call this%dftbp%setupCalculator(input)
364 call this%dftbp%setGeometry(this%coords)
365#endif
366
367 pop_sub(dftb_init)
368 end subroutine dftb_init
369
370 ! ---------------------------------------------------------
371 subroutine dftb_init_interaction(this, interaction)
372 class(dftb_t), target, intent(inout) :: this
373 class(interaction_t), intent(inout) :: interaction
374
375 push_sub(dftb_init_interaction)
376
377 select type (interaction)
378 class default
379 message(1) = "Trying to initialize an unsupported interaction by DFTB+."
380 call messages_fatal(1, namespace=this%namespace)
381 end select
382
383 pop_sub(dftb_init_interaction)
384 end subroutine dftb_init_interaction
385
386 ! ---------------------------------------------------------
387 subroutine dftb_initialize(this)
388 class(dftb_t), intent(inout) :: this
389
390 push_sub(dftb_initialize)
391
392#ifdef HAVE_DFTBPLUS
393 select type (algo => this%algo)
394 class is (propagator_t)
395 select case (this%dynamics)
396 case (bo)
397 call this%dftbp%getGradients(this%gradients)
398 this%tot_force = -this%gradients
399 case (ehrenfest)
400 call this%dftbp%getEnergy(this%energy)
401 call this%dftbp%initializeTimeProp(algo%dt, this%laser_field, .false.)
402 end select
403 end select
404#endif
405
406 pop_sub(dftb_initialize)
407 end subroutine dftb_initialize
408
409 ! ---------------------------------------------------------
410 logical function dftb_do_algorithmic_operation(this, operation, updated_quantities) result(done)
411 class(dftb_t), intent(inout) :: this
412 class(algorithmic_operation_t), intent(in) :: operation
413 character(len=:), allocatable, intent(out) :: updated_quantities(:)
414
415 integer :: ii, jj, il
416 type(tdf_t) :: ff, phi
417 complex(real64) :: amp, pol(3)
418 real(real64) :: time, omega
419
421 call profiling_in(trim(this%namespace%get())//":"//trim(operation%id))
422
423 select type (algo => this%algo)
424 class is (propagator_t)
425
426 done = .true.
427 select case (this%dynamics)
428 case (bo)
429 ! Born-Oppenheimer dynamics
430 select case (operation%id)
432 ! Do nothing
433
434 case (verlet_start)
435 safe_allocate(this%prev_acc(1:3, 1:this%n_atom, 1))
436 do jj = 1, this%n_atom
437 this%acc(1:3, jj) = this%tot_force(1:3, jj) / this%mass(jj)
438 end do
439
440 case (verlet_finish)
441 safe_deallocate_a(this%prev_acc)
442
443 case (verlet_update_pos)
444 do jj = 1, this%n_atom
445 this%coords(1:3, jj) = this%coords(1:3, jj) + algo%dt * this%vel(1:3, jj) &
446 + m_half * algo%dt**2 * this%acc(1:3, jj)
447 end do
448 updated_quantities = ["position"]
449
450 case (verlet_compute_acc)
451 do ii = size(this%prev_acc, dim=3) - 1, 1, -1
452 this%prev_acc(1:3, 1:this%n_atom, ii + 1) = this%prev_acc(1:3, 1:this%n_atom, ii)
453 end do
454 this%prev_acc(1:3, 1:this%n_atom, 1) = this%acc(1:3, 1:this%n_atom)
455#ifdef HAVE_DFTBPLUS
456 call this%dftbp%setGeometry(this%coords)
457 call this%dftbp%getGradients(this%gradients)
458 this%tot_force = -this%gradients
459#endif
460 do jj = 1, this%n_atom
461 this%acc(1:3, jj) = this%tot_force(1:3, jj) / this%mass(jj)
462 end do
463
465 this%vel(1:3, 1:this%n_atom) = this%vel(1:3, 1:this%n_atom) &
466 + m_half * algo%dt * (this%prev_acc(1:3, 1:this%n_atom, 1) + &
467 this%acc(1:3, 1:this%n_atom))
468 updated_quantities = ["velocity"]
469
470 case default
471 done = .false.
472 end select
473
474 case (ehrenfest)
475 ! Ehrenfest dynamics
476 select case (operation%id)
478 ! Do nothing
479 case (verlet_start)
480 !Do nothing
481 case (verlet_finish)
482 !Do nothing
483 case (verlet_update_pos)
484 this%field = m_zero
485 time = this%iteration%value()
486 do il = 1, this%lasers%no_lasers
487 ! get properties of laser
488 call laser_get_f(this%lasers%lasers(il), ff)
489 call laser_get_phi(this%lasers%lasers(il), phi)
490 omega = laser_carrier_frequency(this%lasers%lasers(il))
491 pol = laser_polarization(this%lasers%lasers(il))
492 ! calculate electric field from laser
493 amp = tdf(ff, time) * exp(m_zi * (omega*time + tdf(phi, time)))
494 this%field(1:3) = this%field(1:3) + real(amp*pol(1:3), real64)
495 end do
496#ifdef HAVE_DFTBPLUS
497 call this%dftbp%setTdElectricField(this%field)
498 call this%dftbp%doOneTdStep(this%iteration%counter(), atomnetcharges=this%atom_charges, coord=this%coords,&
499 force=this%tot_force, energy=this%energy)
500#endif
501 case (verlet_compute_acc)
502 !Do nothing
504 !Do nothing
505 case default
506 done = .false.
507 end select
508
509 end select
510
511 end select
512
513 call profiling_out(trim(this%namespace%get())//":"//trim(operation%id))
516
517 ! ---------------------------------------------------------
518 logical function dftb_is_tolerance_reached(this, tol) result(converged)
519 class(dftb_t), intent(in) :: this
520 real(real64), intent(in) :: tol
521
523
524 ! this routine is never called at present, no reason to be here
525 assert(.false.)
526 converged = .false.
527
529 end function dftb_is_tolerance_reached
530
531 ! ---------------------------------------------------------
532 subroutine dftb_output_start(this)
533 class(dftb_t), intent(inout) :: this
534
535 push_sub(dftb_output_start)
536
537 select type (algo => this%algo)
538 class is (propagator_t)
539 ! Create output handle
540 call io_mkdir('td.general', this%namespace)
541 if (mpi_grp_is_root(mpi_world)) then
542 call write_iter_init(this%output_handle(output_coordinates), 0, algo%dt, &
543 trim(io_workpath("td.general/coordinates", this%namespace)))
544 call write_iter_init(this%output_handle(output_forces), 0, algo%dt, &
545 trim(io_workpath("td.general/forces", this%namespace)))
546 end if
547
548 ! Output info for first iteration
549 call this%output_write()
550 end select
551
552 pop_sub(dftb_output_start)
553 end subroutine dftb_output_start
554
555 ! ---------------------------------------------------------
556 subroutine dftb_output_finish(this)
557 class(dftb_t), intent(inout) :: this
558
559 push_sub(dftb_output_finish)
560
561 select type (algo => this%algo)
562 class is (propagator_t)
563 if (mpi_grp_is_root(mpi_world)) then
564 call write_iter_end(this%output_handle(output_coordinates))
565 call write_iter_end(this%output_handle(output_forces))
566 end if
567 end select
568
569 pop_sub(dftb_output_finish)
570 end subroutine dftb_output_finish
571
572 ! ---------------------------------------------------------
573 subroutine dftb_output_write(this)
574 class(dftb_t), intent(inout) :: this
575
576 integer :: idir, iat, iout
577 character(len=50) :: aux
578 character(1) :: out_label(2)
579 real(real64) :: tmp(3)
580
581 if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
582
583 push_sub(dftb_output_write)
584 call profiling_in(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
585
586 select type (algo => this%algo)
587 class is (propagator_t)
588 out_label(1) = "x"
589 out_label(2) = "f"
590
591 if (this%iteration%counter() == 0) then
592 ! header
593 do iout = 1, 2
594 call write_iter_clear(this%output_handle(iout))
595 call write_iter_string(this%output_handle(iout),'#####################################################################')
596 call write_iter_nl(this%output_handle(iout))
597 call write_iter_string(this%output_handle(iout),'# HEADER')
598 call write_iter_nl(this%output_handle(iout))
599
600 ! first line: column names
601 call write_iter_header_start(this%output_handle(iout))
602
603 do iat = 1, this%n_atom
604 do idir = 1, 3
605 write(aux, '(a1,a1,i3,a1,i3,a1)') out_label(iout),'(', iat, ',', idir, ')'
606 call write_iter_header(this%output_handle(iout), aux)
607 end do
608 end do
609 call write_iter_nl(this%output_handle(iout))
610
611 ! second line: units
612 call write_iter_string(this%output_handle(iout), '#[Iter n.]')
613 call write_iter_header(this%output_handle(iout), '[' // trim(units_abbrev(units_out%time)) // ']')
614 end do
615
616 do iat = 1, this%n_atom
617 do idir = 1, 3
618 call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%length)) // ']')
619 call write_iter_header(this%output_handle(output_forces), '[' // trim(units_abbrev(units_out%force)) // ']')
620 end do
621 end do
622
623 do iout = 1, 2
624 call write_iter_nl(this%output_handle(iout))
625 call write_iter_string(this%output_handle(iout),'#######################################################################')
626 call write_iter_nl(this%output_handle(iout))
627 end do
628 end if
629
630 call write_iter_start(this%output_handle(output_coordinates))
631 call write_iter_start(this%output_handle(output_forces))
632
633 do iat = 1, this%n_atom
634 ! Position
635 tmp(1:3) = units_from_atomic(units_out%length, this%coords(1:3, iat))
636 call write_iter_double(this%output_handle(output_coordinates), tmp, 3)
637 ! Force
638 tmp(1:3) = units_from_atomic(units_out%force, this%tot_force(1:3, iat))
639 call write_iter_double(this%output_handle(output_forces), tmp, 3)
640 end do
641
642 call write_iter_nl(this%output_handle(output_coordinates))
643 call write_iter_nl(this%output_handle(output_forces))
644
645 end select
646
647 call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
648 pop_sub(dftb_output_write)
649 end subroutine dftb_output_write
650
651 ! ---------------------------------------------------------
652 subroutine dftb_init_interaction_as_partner(partner, interaction)
653 class(dftb_t), intent(in) :: partner
654 class(interaction_surrogate_t), intent(inout) :: interaction
655
657
658 select type (interaction)
659 class default
660 message(1) = "Unsupported interaction."
661 call messages_fatal(1, namespace=partner%namespace)
662 end select
663
667 ! ---------------------------------------------------------
668 subroutine dftb_copy_quantities_to_interaction(partner, interaction)
669 class(dftb_t), intent(inout) :: partner
670 class(interaction_surrogate_t), intent(inout) :: interaction
671
673
674 select type (interaction)
675 class default
676 message(1) = "Unsupported interaction."
677 call messages_fatal(1, namespace=partner%namespace)
678 end select
679
682
683 ! ---------------------------------------------------------
684 subroutine dftb_update_interactions_start(this)
685 class(dftb_t), intent(inout) :: this
686
688
689 ! Store previous force, as it is used as SCF criterium
690 this%prev_tot_force(1:3, 1:this%n_atom) = this%tot_force(1:3, 1:this%n_atom)
691
693 end subroutine dftb_update_interactions_start
694
695 ! ---------------------------------------------------------
696 subroutine dftb_update_kinetic_energy(this)
697 class(dftb_t), intent(inout) :: this
698
700
701 this%kinetic_energy = m_zero
702
704 end subroutine dftb_update_kinetic_energy
705
706 ! ---------------------------------------------------------
707 subroutine dftb_restart_write_data(this)
708 class(dftb_t), intent(inout) :: this
709
711
712 message(1) = "DFTB system "//trim(this%namespace%get())//" cannot write restart data."
713 call messages_warning(1, namespace=this%namespace)
714
716 end subroutine dftb_restart_write_data
717
718 ! ---------------------------------------------------------
719 ! this function returns true if restart data could be read
720 logical function dftb_restart_read_data(this)
721 class(dftb_t), intent(inout) :: this
722
723 push_sub(dftb_restart_read_data)
724
725 ! restarting not yet supported
726 dftb_restart_read_data = .false.
727
729 end function dftb_restart_read_data
730
731 ! ---------------------------------------------------------
732 subroutine dftb_finalize(this)
733 type(dftb_t), intent(inout) :: this
734
735 push_sub(dftb_finalize)
736
737 safe_deallocate_a(this%coords)
738 safe_deallocate_a(this%acc)
739 safe_deallocate_a(this%vel)
740 safe_deallocate_a(this%tot_force)
741 safe_deallocate_a(this%prev_tot_force)
742 safe_deallocate_a(this%gradients)
743 safe_deallocate_a(this%species)
744 safe_deallocate_a(this%mass)
745 call ion_dynamics_end(this%ions_dyn)
746
747 deallocate(this%ions)
748
749 ! No call to safe_deallocate macro here, as it gives an ICE with gfortran
750 if (associated(this%lasers)) then
751 deallocate(this%lasers)
752 end if
753
754#ifdef HAVE_DFTBPLUS
755 call tdftbplus_destruct(this%dftbp)
756#endif
757
758 call system_end(this)
759
760 pop_sub(dftb_finalize)
761 end subroutine dftb_finalize
762
763end module dftb_oct_m
764
765!! Local Variables:
766!! mode: f90
767!! coding: utf-8
768!! End:
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:180
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
Definition: write_iter.F90:167
double exp(double __x) __attribute__((__nothrow__
This module implements the basic elements defining algorithms.
Definition: algorithm.F90:141
subroutine dftb_output_write(this)
Definition: dftb.F90:667
subroutine dftb_update_interactions_start(this)
Definition: dftb.F90:778
subroutine dftb_init_interaction_as_partner(partner, interaction)
Definition: dftb.F90:746
subroutine dftb_output_finish(this)
Definition: dftb.F90:650
subroutine dftb_restart_write_data(this)
Definition: dftb.F90:801
subroutine dftb_init_interaction(this, interaction)
Definition: dftb.F90:465
logical function dftb_restart_read_data(this)
Definition: dftb.F90:814
logical function dftb_do_algorithmic_operation(this, operation, updated_quantities)
Definition: dftb.F90:504
class(dftb_t) function, pointer dftb_constructor(namespace)
The factory routine (or constructor) allocates a pointer of the corresponding type and then calls the...
Definition: dftb.F90:216
subroutine, public dftb_init(this, namespace)
The init routine is a module level procedure This has the advantage that different classes can have d...
Definition: dftb.F90:240
subroutine dftb_update_kinetic_energy(this)
Definition: dftb.F90:790
subroutine dftb_initialize(this)
Definition: dftb.F90:481
integer, parameter output_forces
Definition: dftb.F90:152
subroutine dftb_finalize(this)
Definition: dftb.F90:826
logical function dftb_is_tolerance_reached(this, tol)
Definition: dftb.F90:612
subroutine dftb_output_start(this)
Definition: dftb.F90:626
integer, parameter bo
Definition: dftb.F90:204
subroutine dftb_copy_quantities_to_interaction(partner, interaction)
Definition: dftb.F90:762
real(real64), parameter, public m_zero
Definition: global.F90:188
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
This module defines the abstract interaction_t class, and some auxiliary classes for interactions.
Definition: io.F90:114
character(len=max_path_len) function, public io_workpath(path, namespace)
Definition: io.F90:270
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:311
logical pure function, public ion_dynamics_ions_move(this)
subroutine, public ion_dynamics_unfreeze(this)
Unfreezes the ionic movement.
subroutine, public ion_dynamics_init(this, namespace, ions)
subroutine, public ion_dynamics_end(this)
complex(real64) function, dimension(3), public laser_polarization(laser)
Definition: lasers.F90:726
subroutine, public lasers_parse_external_fields(this)
Definition: lasers.F90:244
subroutine, public laser_get_f(laser, ff)
Definition: lasers.F90:767
real(real64) function, public laser_carrier_frequency(laser)
Definition: lasers.F90:621
subroutine, public laser_get_phi(laser, phi)
Definition: lasers.F90:795
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
character(len=512), private msg
Definition: messages.F90:165
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
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_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:430
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:266
This module implements the multisystem debug functionality.
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
This module implements the basic propagator framework.
Definition: propagator.F90:117
character(len=algo_label_len), parameter, public store_current_status
Definition: propagator.F90:169
character(len=30), parameter, public verlet_start
character(len=30), parameter, public verlet_compute_acc
character(len=30), parameter, public verlet_update_pos
character(len=30), parameter, public verlet_finish
character(len=30), parameter, public verlet_compute_vel
This module defines the quantity_t class and the IDs for quantities, which can be exposed by a system...
Definition: quantity.F90:138
This module implements the abstract system type.
Definition: system.F90:118
subroutine, public system_end(this)
Definition: system.F90:1145
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.
Explicit interfaces to C functions, defined in write_iter_low.cc.
Definition: write_iter.F90:114
Descriptor of one algorithmic operation.
Definition: algorithm.F90:163
class for a tight binding
Definition: dftb.F90:158
abstract interaction class
surrogate interaction class to avoid circular dependencies between modules.
Abstract class implementing propagators.
Definition: propagator.F90:138
Systems (system_t) can expose quantities that can be used to calculate interactions with other system...
Definition: quantity.F90:171
Abstract class for systems.
Definition: system.F90:172
int true(void)
subroutine input()
Definition: vdw.F90:255