33 use,
intrinsic :: iso_fortran_env
68 integer,
parameter :: &
69 OUTPUT_COORDINATES = 1, &
74 type(c_ptr),
public :: output_handle(2)
75 type(space_t),
public :: space
77 real(real64),
allocatable,
public :: mass(:)
78 real(real64),
allocatable,
public :: pos(:,:)
79 real(real64),
allocatable,
public :: vel(:,:)
80 real(real64),
allocatable,
public :: tot_force(:,:)
81 real(real64),
allocatable,
public :: lj_epsilon(:)
82 real(real64),
allocatable,
public :: lj_sigma(:)
83 logical,
allocatable,
public :: fixed(:)
84 type(propagator_data_t),
public :: prop_data
114 class(classical_particles_t),
intent(inout) :: this
115 integer,
intent(in) :: np
120 safe_allocate(this%mass(1:np))
122 safe_allocate(this%pos(1:this%space%dim, 1:np))
124 safe_allocate(this%vel(1:this%space%dim, 1:np))
126 safe_allocate(this%tot_force(1:this%space%dim, 1:np))
128 safe_allocate(this%fixed(1:np))
130 safe_allocate(this%lj_epsilon(1:np))
132 safe_allocate(this%lj_sigma(1:np))
137 allocate(this%supported_interactions(0))
138 allocate(this%supported_interactions_as_partner(0))
141 call this%quantities%add(
quantity_t(
"position", updated_on_demand = .false.))
142 call this%quantities%add(
quantity_t(
"velocity", updated_on_demand = .false.))
145 call this%quantities%add(
quantity_t(
"mass", updated_on_demand = .false., always_available=.
true.))
152 class(classical_particles_t),
intent(out) :: this
153 class(classical_particles_t),
intent(in) :: cp_in
158 safe_allocate_source_a(this%mass, cp_in%mass)
159 safe_allocate_source_a(this%pos, cp_in%pos)
160 safe_allocate_source_a(this%vel, cp_in%vel)
161 safe_allocate_source_a(this%tot_force, cp_in%tot_force)
162 safe_allocate_source_a(this%fixed, cp_in%fixed)
164 this%kinetic_energy = cp_in%kinetic_energy
166 this%quantities = cp_in%quantities
167 this%supported_interactions = cp_in%supported_interactions
169 this%prop_data = cp_in%prop_data
181 select type (interaction)
183 message(1) =
"Trying to initialize an unsupported interaction by classical particles."
194 character(len=:),
allocatable,
intent(out) :: updated_quantities(:)
197 real(real64),
allocatable :: tmp_pos(:,:,:), tmp_vel(:,:,:)
198 real(real64) :: factor
201 call profiling_in(trim(this%namespace%get())//
":"//trim(operation%id))
204 select type (prop => this%algo)
207 select case (operation%id)
209 this%prop_data%save_pos(:, 1:this%np) = this%pos(:, 1:this%np)
210 this%prop_data%save_vel(:, 1:this%np) = this%vel(:, 1:this%np)
213 if (.not. this%prop_data%initialized)
then
214 call this%prop_data%initialize(prop, this%space%dim, this%np)
216 if (this%fixed(ip))
then
217 this%prop_data%acc(:, ip) =
m_zero
219 this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
225 call this%prop_data%end()
228 call this%prop_data%end()
231 this%pos(:, 1:this%np) = this%pos(:, 1:this%np) + prop%dt * this%vel(:, 1:this%np) &
232 +
m_half * prop%dt**2 * this%prop_data%acc(:, 1:this%np)
233 updated_quantities = [
"position"]
236 do ii =
size(this%prop_data%prev_acc, dim=3) - 1, 1, -1
237 this%prop_data%prev_acc(:, 1:this%np, ii + 1) = this%prop_data%prev_acc(:, 1:this%np, ii)
240 this%prop_data%prev_acc(:, ip, 1) = this%prop_data%acc(:, ip)
241 if (this%fixed(ip))
then
242 this%prop_data%acc(:, ip) =
m_zero
244 this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
249 this%vel(:, 1:this%np) = this%vel(:, 1:this%np) &
250 +
m_half * prop%dt * (this%prop_data%prev_acc(:, 1:this%np, 1) + this%prop_data%acc(:, 1:this%np))
251 updated_quantities = [
"velocity"]
254 if (.not. this%prop_data%initialized)
then
255 call this%prop_data%initialize(prop, this%space%dim, this%np)
257 if (this%fixed(ip))
then
258 this%prop_data%acc(:, ip) =
m_zero
260 this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
262 this%prop_data%prev_acc(:, ip, 1) = this%prop_data%acc(:, ip)
267 this%pos(:, 1:this%np) = this%pos(:, 1:this%np) + prop%dt * this%vel(:, 1:this%np) + &
268 m_one/6.0_real64 * prop%dt**2 * (
m_four*this%prop_data%acc(:, 1:this%np) - this%prop_data%prev_acc(:, 1:this%np, 1))
270 if (.not. prop%predictor_corrector)
then
271 updated_quantities = [
"position"]
275 this%vel(:, 1:this%np) = this%vel(:, 1:this%np) + &
276 m_one/6.0_real64 * prop%dt * (
m_two * this%prop_data%acc(:, 1:this%np) + &
277 5.0_real64 * this%prop_data%prev_acc(:, 1:this%np, 1) - this%prop_data%prev_acc(:, 1:this%np, 2))
278 updated_quantities = [
"velocity"]
281 this%pos(:, 1:this%np) = this%prop_data%save_pos(:, 1:this%np) &
282 + prop%dt * this%prop_data%save_vel(:, 1:this%np) &
283 +
m_one/6.0_real64 * prop%dt**2 * (this%prop_data%acc(:, 1:this%np) &
284 +
m_two * this%prop_data%prev_acc(:, 1:this%np, 1))
285 updated_quantities = [
"position"]
288 this%vel(:, 1:this%np) = this%prop_data%save_vel(:, 1:this%np) &
289 +
m_half * prop%dt * (this%prop_data%acc(:, 1:this%np) + this%prop_data%prev_acc(:, 1:this%np, 1))
290 updated_quantities = [
"velocity"]
293 if (.not. this%prop_data%initialized)
then
294 call this%prop_data%initialize(prop, this%space%dim, this%np)
295 this%prop_data%prev_pos(:, 1:this%np, 1) = this%pos(:, 1:this%np)
296 this%prop_data%prev_vel(:, 1:this%np, 1) = this%vel(:, 1:this%np)
300 call this%prop_data%end()
303 this%pos(:, 1:this%np) = 1.5_real64*this%prop_data%save_pos(:, 1:this%np) &
304 -
m_half*this%prop_data%prev_pos(:, 1:this%np, 1)
305 this%vel(:, 1:this%np) = 1.5_real64*this%prop_data%save_vel(:, 1:this%np) &
306 -
m_half*this%prop_data%prev_vel(:, 1:this%np, 1)
307 this%prop_data%prev_pos(:, 1:this%np, 1) = this%prop_data%save_pos(:, 1:this%np)
308 this%prop_data%prev_vel(:, 1:this%np, 1) = this%prop_data%save_vel(:, 1:this%np)
309 updated_quantities = [
"position",
"velocity"]
313 if (this%fixed(ip))
then
314 this%prop_data%hamiltonian_elements(:, ip) =
m_zero
316 this%prop_data%hamiltonian_elements(:, ip) = this%tot_force(:, ip) / (this%mass(ip) * this%pos(:, ip) +
m_epsilon)
321 safe_allocate(tmp_pos(1:this%space%dim, 1:this%np, 1:2))
322 safe_allocate(tmp_vel(1:this%space%dim, 1:this%np, 1:2))
328 tmp_pos(:, 1:this%np, 1) = this%prop_data%save_pos(:, 1:this%np)
329 tmp_vel(:, 1:this%np, 1) = this%prop_data%save_vel(:, 1:this%np)
330 this%pos(:, 1:this%np) = this%prop_data%save_pos(:, 1:this%np)
331 this%vel(:, 1:this%np) = this%prop_data%save_vel(:, 1:this%np)
335 factor = factor * prop%dt / ii
338 tmp_pos(:, ip, 2) = tmp_vel(:, ip, 1)
339 tmp_vel(:, ip, 2) = this%prop_data%hamiltonian_elements(:, ip) * tmp_pos(:, ip, 1)
341 tmp_pos(:, ip, 1) = tmp_pos(:, ip, 2)
342 tmp_vel(:, ip, 1) = tmp_vel(:, ip, 2)
344 this%pos(:, ip) = this%pos(:, ip) + factor * tmp_pos(:, ip, 1)
345 this%vel(:, ip) = this%vel(:, ip) + factor * tmp_vel(:, ip, 1)
348 safe_deallocate_a(tmp_pos)
349 safe_deallocate_a(tmp_vel)
350 updated_quantities = [
"position",
"velocity"]
354 if (.not. this%is_tolerance_reached(prop%scf_tol))
then
355 this%pos(:, 1:this%np) =
m_half*(this%pos(:, 1:this%np) + this%prop_data%save_pos(:, 1:this%np))
356 this%vel(:, 1:this%np) =
m_half*(this%vel(:, 1:this%np) + this%prop_data%save_vel(:, 1:this%np))
357 updated_quantities = [
"position",
"velocity"]
368 call profiling_out(trim(this%namespace%get())//
":"//trim(operation%id))
375 real(real64),
intent(in) :: tol
378 real(real64) :: change, max_change
385 change = sum((this%prop_data%prev_tot_force(1:this%space%dim, ip) - this%tot_force(1:this%space%dim, ip))**2) / &
387 if (change > max_change)
then
391 converged = max_change < tol**2
393 write(
message(1),
'(a, e13.6, a, e13.6)')
"Debug: -- Maximum change in acceleration ", &
394 sqrt(max_change),
" and tolerance ", tol
403 character(len=*),
intent(in) :: label
409 message(1) =
"Incompatible quantity."
423 select type (interaction)
425 message(1) =
"Unsupported interaction."
439 select type (interaction)
441 message(1) =
"Unsupported interaction."
455 select type (prop => this%algo)
457 if (prop%predictor_corrector)
then
458 this%prop_data%prev_tot_force(1:this%space%dim, 1:this%np) = this%tot_force(1:this%space%dim, 1:this%np)
474 this%tot_force(1:this%space%dim, 1:this%np) =
m_zero
475 call iter%start(this%interactions)
476 do while (iter%has_next())
477 select type (interaction => iter%get_next())
479 this%tot_force(1:this%space%dim, 1:this%np) = this%tot_force(1:this%space%dim, 1:this%np) + &
480 interaction%force(1:this%space%dim, 1:this%np)
496 select type (algo => this%algo)
501 iteration = this%iteration%counter()
502 if (iteration > 0)
then
504 iteration = iteration + 1
507 call io_mkdir(
'td.general', this%namespace)
509 call write_iter_init(this%output_handle(output_coordinates), iteration, dt, &
512 trim(
io_workpath(
"td.general/energy", this%namespace)))
516 if (iteration == 0)
then
517 call this%output_write()
541 integer :: idir, ii, ip
542 character(len=50) :: aux
543 real(real64) :: tmp(this%space%dim)
549 if (this%iteration%counter() == 0)
then
553 '############################################################################')
561 do idir = 1, this%space%dim
562 write(aux,
'(a2,i3,a1,i3,a1)')
'x(', ip,
',', idir,
')'
567 do idir = 1, this%space%dim
568 write(aux,
'(a2,i3,a1,i3,a1)')
'v(', ip,
',', idir,
')'
573 do idir = 1, this%space%dim
574 write(aux,
'(a2,i3,a1,i3,a1)')
'f(', ip,
',', idir,
')'
584 do idir = 1, this%space%dim
589 do idir = 1, this%space%dim
594 do idir = 1, this%space%dim
601 '############################################################################')
610 call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
615 call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
620 call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
626 if (this%iteration%counter() == 0)
then
630 '############################################################################')
652 '############################################################################')
670 integer :: restart_file_unit
679 call profiling_in(trim(this%namespace%get())//
":"//
"RESTART_WRITE")
682 restart_file_unit =
io_open(
'restart/'//
td_dir//
'restart_classical_particles', this%namespace, action=
'write')
683 write(restart_file_unit, *) this%np
684 write(restart_file_unit, *) this%mass(:)
685 write(restart_file_unit, *) this%pos(:,:)
686 write(restart_file_unit, *) this%vel(:,:)
687 write(restart_file_unit, *) this%tot_force(:,:)
690 if (this%iteration%counter() > 0)
then
692 call this%prop_data%restart_write(this%namespace, this%algo)
695 message(1) =
"Successfully wrote restart data for system "//trim(this%namespace%get())
698 call profiling_out(trim(this%namespace%get())//
":"//
"RESTART_WRITE")
706 integer :: restart_file_unit
709 call profiling_in(trim(this%namespace%get())//
":"//
"RESTART_READ")
712 restart_file_unit =
io_open(
'restart/'//
td_dir//
'restart_classical_particles', this%namespace, action=
'read', die=.false.)
713 if (restart_file_unit /= -1)
then
714 read(restart_file_unit, *) this%np
715 read(restart_file_unit, *) this%mass(:)
716 read(restart_file_unit, *) this%pos(:,:)
717 read(restart_file_unit, *) this%vel(:,:)
718 read(restart_file_unit, *) this%tot_force(:,:)
720 call this%prop_data%initialize(this%algo, this%space%dim, this%np)
721 if (this%iteration%counter() > 0)
then
733 message(1) =
"Successfully read restart data for system "//trim(this%namespace%get())
737 message(1) =
"Failed to read restart data for system "//trim(this%namespace%get())
742 call profiling_out(trim(this%namespace%get())//
":"//
"RESTART_READ")
754 this%kinetic_energy =
m_zero
756 this%kinetic_energy = this%kinetic_energy +
m_half * this%mass(ip) * sum(this%vel(:, ip)**2)
765 logical,
optional,
intent(in) :: mask(:)
766 logical,
optional,
intent(in) :: pseudo
767 real(real64) :: pos(this%space%dim)
769 real(real64) :: mass, total_mass
774 assert(.not. this%space%is_periodic())
780 if (
present(mask))
then
781 if (.not. mask(ip)) cycle
786 pos = pos + mass*this%pos(:, ip)
787 total_mass = total_mass + mass
797 real(real64) :: vel(this%space%dim)
799 real(real64) :: mass, total_mass
808 total_mass = total_mass + mass
809 vel = vel + mass*this%vel(:, ip)
819 real(real64) :: pos(this%space%dim)
821 real(real64) :: xmin(this%space%dim), xmax(this%space%dim)
829 do idir = 1, this%space%dim
830 if (this%pos(idir, ip) > xmax(idir)) xmax(idir) = this%pos(idir, ip)
831 if (this%pos(idir, ip) < xmin(idir)) xmin(idir) = this%pos(idir, ip)
835 pos = (xmax + xmin)/
m_two
843 real(real64),
intent(out) :: x(this%space%dim)
844 real(real64),
intent(out) :: x2(this%space%dim)
847 real(real64) :: rmax, r, r2
851 assert(.not. this%space%is_periodic())
856 do jp = 1, this%np/2 + 1
857 r = norm2(this%pos(:, ip) - this%pos(:, jp))
860 x = this%pos(:, ip) - this%pos(:, jp)
869 r2 = sum(x * this%pos(:, ip))
870 r = norm2(this%pos(:, ip) - r2*x)
873 x2 = this%pos(:, ip) - r2*x
885 real(real64),
intent(out) :: x(this%space%dim)
886 real(real64),
intent(out) :: x2(this%space%dim)
887 logical,
intent(in) :: pseudo
889 real(real64) :: mass, tinertia(this%space%dim, this%space%dim), eigenvalues(this%space%dim)
890 integer :: ii, jj, ip
895 assert(.not. this%space%is_periodic())
901 if (.not. pseudo) mass = this%mass(ip)
902 do ii = 1, this%space%dim
903 tinertia(ii, :) = tinertia(ii, :) - mass*this%pos(ii, ip)*this%pos(:, ip)
904 tinertia(ii, ii) = tinertia(ii, ii) + mass*sum(this%pos(:, ip)**2)
911 write(
message(1),
'(a)')
'Moment of pseudo-inertia tensor [' // trim(
units_abbrev(unit)) //
']'
913 write(
message(1),
'(a)')
'Moment of inertia tensor [amu*' // trim(
units_abbrev(unit)) //
']'
916 call output_tensor(tinertia, this%space%dim, unit, write_average = .
true., namespace=this%namespace)
925 jj = maxloc(abs(tinertia(:,ii)), dim = 1)
926 if (tinertia(jj,ii) <
m_zero) tinertia(:,ii) = -tinertia(:,ii)
940 safe_deallocate_a(this%mass)
941 safe_deallocate_a(this%pos)
942 safe_deallocate_a(this%vel)
943 safe_deallocate_a(this%tot_force)
944 safe_deallocate_a(this%fixed)
945 safe_deallocate_a(this%lj_epsilon)
946 safe_deallocate_a(this%lj_sigma)
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
This module implements the basic elements defining algorithms.
subroutine classical_particles_output_finish(this)
subroutine classical_particles_axis_large(this, x, x2)
logical function classical_particles_is_tolerance_reached(this, tol)
subroutine classical_particles_output_write(this)
subroutine, public classical_particles_update_quantity(this, label)
integer, parameter output_energy
logical function classical_particles_do_algorithmic_operation(this, operation, updated_quantities)
subroutine classical_particles_update_kinetic_energy(this)
subroutine classical_particles_axis_inertia(this, x, x2, pseudo)
This subroutine assumes that the origin of the coordinates is the center of mass of the system.
subroutine classical_particles_output_start(this)
logical function classical_particles_restart_read_data(this)
subroutine, public classical_particles_init_interaction_as_partner(partner, interaction)
subroutine, public classical_particles_end(this)
subroutine, public classical_particles_init(this, np)
The init routine is a module level procedure This has the advantage that different classes can have d...
subroutine, public classical_particles_copy(this, cp_in)
subroutine classical_particles_update_interactions_finish(this)
subroutine classical_particles_update_interactions_start(this)
subroutine, public classical_particles_copy_quantities_to_interaction(partner, interaction)
real(real64) function, dimension(this%space%dim) classical_particles_center_of_mass_vel(this)
subroutine, public classical_particles_init_interaction(this, interaction)
subroutine, public classical_particles_restart_write_data(this)
real(real64) function, dimension(this%space%dim) classical_particles_center_of_mass(this, mask, pseudo)
real(real64) function, dimension(this%space%dim) classical_particles_center(this)
real(real64), parameter, public m_two
real(real64), parameter, public m_huge
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_epsilon
character(len= *), parameter, public td_dir
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.
subroutine, public io_close(iunit, grp)
character(len=max_path_len) function, public io_workpath(path, namespace)
subroutine, public io_mkdir(fname, namespace, parents)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
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_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.
This module implements the multisystem debug functionality.
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.
character(len=algo_label_len), parameter, public beeman_start
character(len=algo_label_len), parameter, public beeman_finish
character(len=algo_label_len), parameter, public beeman_compute_acc
character(len=algo_label_len), parameter, public beeman_correct_pos
character(len=algo_label_len), parameter, public beeman_predict_pos
character(len=algo_label_len), parameter, public beeman_correct_vel
character(len=algo_label_len), parameter, public beeman_predict_vel
character(len=algo_label_len), parameter, public expmid_2step_predict_dt
character(len=algo_label_len), parameter, public expmid_2step_finish
character(len=algo_label_len), parameter, public expmid_2step_correct_dt_2
character(len=algo_label_len), parameter, public update_hamiltonian
character(len=algo_label_len), parameter, public expmid_2step_predict_dt_2
character(len=algo_label_len), parameter, public expmid_2step_start
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
This module is intended to contain simple general-purpose utility functions and procedures.
subroutine, public output_tensor(tensor, ndim, unit, write_average, iunit, namespace)
Explicit interfaces to C functions, defined in write_iter_low.cc.
Descriptor of one algorithmic operation.
These class extend the list and list iterator to make an interaction list.
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.