34 use,
intrinsic :: iso_fortran_env
59 integer,
parameter :: &
60 OUTPUT_COORDINATES = 1, &
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(:)
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
87 type(TDftbPlus) :: dftbp
107 procedure dftb_constructor
111 integer,
parameter :: &
123 class(dftb_t),
pointer :: sys
124 type(namespace_t),
intent(in) :: namespace
129 message(1) =
"DFTB+ system not available. This feature requires compiling Octopus with the DFTB+ library"
147 class(dftb_t),
target,
intent(inout) :: this
148 type(namespace_t),
intent(in) :: namespace
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
156 real(real64) :: initial_temp
160 type(fnode),
pointer :: proot, pgeo, pham, pdftb, pmaxang, pslakos, ptype2files, panalysis
161 type(fnode),
pointer :: pparseropts
162 type(fnode),
pointer :: pelecdyn, pperturb, plaser
167 this%namespace = namespace
170 allocate(this%supported_interactions(0))
171 allocate(this%supported_interactions_as_partner(0))
174 call this%quantities%add(
quantity_t(
"position", updated_on_demand = .false.))
175 call this%quantities%add(
quantity_t(
"velocity", updated_on_demand = .false.))
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))
195 this%labels(1) = trim(this%ions%atom(1)%label)
197 this%coords = this%ions%pos
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
203 this%labels(ispec) = trim(this%ions%atom(ii)%label)
206 if (trim(this%ions%atom(ii)%label) == trim(this%labels(jj)))
then
207 this%species(ii) = jj
229 if (
parse_block(namespace,
'MaxAngularMomentum', blk) == 0)
then
231 if (n_maxang_block /= this%ions%nspecies)
then
235 do ii = 1, n_maxang_block
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))
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))
266 if (this%dynamics ==
bo)
then
284 if (this%lasers%no_lasers > 0)
then
285 this%laser_field = .
true.
287 this%laser_field = .false.
291 call tdftbplus_init(this%dftbp)
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.)
310 call parse_variable(namespace,
'SccTolerance', 1e-9_real64, this%scc_tolerance)
312 call setchildvalue(pdftb,
"SccTolerance", this%scc_tolerance)
315 call setchild(pdftb,
"MaxAngularMomentum", pmaxang)
317 do ii = 1, this%ions%nspecies
318 call setchildvalue(pmaxang, this%labels(ii), max_ang_mom(ii))
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")
331 call setchild(proot,
"Analysis", panalysis)
332 call setchildvalue(panalysis,
"CalculateForces", .
true.)
334 call setchild(proot,
"ParserOptions", pparseropts)
335 call setchildvalue(pparseropts,
"ParserVersion", 5)
337 if (this%dynamics == ehrenfest)
then
338 call setchild(proot,
"ElectronDynamics", pelecdyn)
341 call setchildvalue(pelecdyn,
"InitialTemperature", initial_temp)
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)
351 call setchildvalue(plaser,
"LaserEnergy",
m_one)
352 call setchildvalue(pelecdyn,
"FieldStrength",
m_one)
354 call setchild(pperturb,
"None", plaser)
358 message(1) =
'Input tree in HSD format:'
360 call dumphsd(
input%hsdTree, stdout)
363 call this%dftbp%setupCalculator(
input)
364 call this%dftbp%setGeometry(this%coords)
372 class(
dftb_t),
target,
intent(inout) :: this
377 select type (interaction)
379 message(1) =
"Trying to initialize an unsupported interaction by DFTB+."
388 class(
dftb_t),
intent(inout) :: this
393 select type (algo => this%algo)
395 select case (this%dynamics)
397 call this%dftbp%getGradients(this%gradients)
398 this%tot_force = -this%gradients
400 call this%dftbp%getEnergy(this%energy)
401 call this%dftbp%initializeTimeProp(algo%dt, this%laser_field, .false.)
411 class(
dftb_t),
intent(inout) :: this
413 character(len=:),
allocatable,
intent(out) :: updated_quantities(:)
415 integer :: ii, jj, il
416 type(
tdf_t) :: ff, phi
417 complex(real64) :: amp, pol(3)
418 real(real64) :: time, omega
421 call profiling_in(trim(this%namespace%get())//
":"//trim(operation%id))
423 select type (algo => this%algo)
427 select case (this%dynamics)
430 select case (operation%id)
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)
441 safe_deallocate_a(this%prev_acc)
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)
448 updated_quantities = [
"position"]
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)
454 this%prev_acc(1:3, 1:this%n_atom, 1) = this%acc(1:3, 1:this%n_atom)
456 call this%dftbp%setGeometry(this%coords)
457 call this%dftbp%getGradients(this%gradients)
458 this%tot_force = -this%gradients
460 do jj = 1, this%n_atom
461 this%acc(1:3, jj) = this%tot_force(1:3, jj) / this%mass(jj)
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"]
476 select case (operation%id)
485 time = this%iteration%value()
486 do il = 1, this%lasers%no_lasers
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)
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)
513 call profiling_out(trim(this%namespace%get())//
":"//trim(operation%id))
519 class(
dftb_t),
intent(in) :: this
520 real(real64),
intent(in) :: tol
533 class(
dftb_t),
intent(inout) :: this
537 select type (algo => this%algo)
540 call io_mkdir(
'td.general', this%namespace)
542 call write_iter_init(this%output_handle(output_coordinates), 0, algo%dt, &
543 trim(
io_workpath(
"td.general/coordinates", this%namespace)))
545 trim(
io_workpath(
"td.general/forces", this%namespace)))
549 call this%output_write()
557 class(
dftb_t),
intent(inout) :: this
561 select type (algo => this%algo)
574 class(
dftb_t),
intent(inout) :: this
576 integer :: idir, iat, iout
577 character(len=50) :: aux
578 character(1) :: out_label(2)
579 real(real64) :: tmp(3)
584 call profiling_in(trim(this%namespace%get())//
":"//
"OUTPUT_WRITE")
586 select type (algo => this%algo)
591 if (this%iteration%counter() == 0)
then
595 call write_iter_string(this%output_handle(iout),
'#####################################################################')
603 do iat = 1, this%n_atom
605 write(aux,
'(a1,a1,i3,a1,i3,a1)') out_label(iout),
'(', iat,
',', idir,
')'
616 do iat = 1, this%n_atom
625 call write_iter_string(this%output_handle(iout),
'#######################################################################')
633 do iat = 1, this%n_atom
647 call profiling_out(trim(this%namespace%get())//
":"//
"OUTPUT_WRITE")
653 class(
dftb_t),
intent(in) :: partner
658 select type (interaction)
660 message(1) =
"Unsupported interaction."
669 class(
dftb_t),
intent(inout) :: partner
674 select type (interaction)
676 message(1) =
"Unsupported interaction."
685 class(
dftb_t),
intent(inout) :: this
690 this%prev_tot_force(1:3, 1:this%n_atom) = this%tot_force(1:3, 1:this%n_atom)
697 class(
dftb_t),
intent(inout) :: this
701 this%kinetic_energy =
m_zero
708 class(
dftb_t),
intent(inout) :: this
712 message(1) =
"DFTB system "//trim(this%namespace%get())//
" cannot write restart data."
721 class(
dftb_t),
intent(inout) :: this
733 type(
dftb_t),
intent(inout) :: this
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)
747 deallocate(this%ions)
750 if (
associated(this%lasers))
then
751 deallocate(this%lasers)
755 call tdftbplus_destruct(this%dftbp)
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
double exp(double __x) __attribute__((__nothrow__
This module implements the basic elements defining algorithms.
subroutine dftb_output_write(this)
subroutine dftb_update_interactions_start(this)
subroutine dftb_init_interaction_as_partner(partner, interaction)
subroutine dftb_output_finish(this)
subroutine dftb_restart_write_data(this)
subroutine dftb_init_interaction(this, interaction)
logical function dftb_restart_read_data(this)
logical function dftb_do_algorithmic_operation(this, operation, updated_quantities)
class(dftb_t) function, pointer dftb_constructor(namespace)
The factory routine (or constructor) allocates a pointer of the corresponding type and then calls the...
subroutine, public dftb_init(this, namespace)
The init routine is a module level procedure This has the advantage that different classes can have d...
subroutine dftb_update_kinetic_energy(this)
subroutine dftb_initialize(this)
integer, parameter output_forces
subroutine dftb_finalize(this)
logical function dftb_is_tolerance_reached(this, tol)
subroutine dftb_output_start(this)
subroutine dftb_copy_quantities_to_interaction(partner, interaction)
real(real64), parameter, public m_zero
complex(real64), parameter, public m_zi
real(real64), parameter, public m_half
real(real64), parameter, public m_one
This module defines the abstract interaction_t class, and some auxiliary classes for interactions.
character(len=max_path_len) function, public io_workpath(path, namespace)
subroutine, public io_mkdir(fname, namespace, parents)
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)
subroutine, public lasers_parse_external_fields(this)
subroutine, public laser_get_f(laser, ff)
real(real64) function, public laser_carrier_frequency(laser)
subroutine, public laser_get_phi(laser, phi)
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
character(len=512), private msg
subroutine, public messages_warning(no_lines, all_nodes, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
type(mpi_grp_t), public mpi_world
This module implements the multisystem debug functionality.
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
This module implements the basic propagator framework.
character(len=algo_label_len), parameter, public store_current_status
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...
This module implements the abstract system type.
subroutine, public system_end(this)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
character(len=20) pure function, public units_abbrev(this)
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.
Descriptor of one algorithmic operation.
class for a tight binding
abstract interaction class
surrogate interaction class to avoid circular dependencies between modules.
Abstract class implementing propagators.
Systems (system_t) can expose quantities that can be used to calculate interactions with other system...
Abstract class for systems.