Octopus
classical_particles.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2019 N. Tancogne-Dejean
3!! Copyright (C) 2020 M. Oliveira
4!! Copyright (C) 2021 S. Ohlmann, I. Albar, A. Obzhirov, A. Geondzhian, M. Lawan
5!! Copyright (C) 2022 M. Oliveira
6!!
7!! This program is free software; you can redistribute it and/or modify
8!! it under the terms of the GNU General Public License as published by
9!! the Free Software Foundation; either version 2, or (at your option)
10!! any later version.
11!!
12!! This program is distributed in the hope that it will be useful,
13!! but WITHOUT ANY WARRANTY; without even the implied warranty of
14!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15!! GNU General Public License for more details.
16!!
17!! You should have received a copy of the GNU General Public License
18!! along with this program; if not, write to the Free Software
19!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20!! 02110-1301, USA.
21!!
22
23#include "global.h"
24
27 use debug_oct_m
29 use global_oct_m
32 use io_oct_m
33 use, intrinsic :: iso_fortran_env
36 use mpi_oct_m
39 use parser_oct_m
47 use space_oct_m
48 use system_oct_m
49 use unit_oct_m
51 use utils_oct_m
53
54 implicit none
55
56 private
57 public :: &
67
68 integer, parameter :: &
69 OUTPUT_COORDINATES = 1, &
71
72 type, extends(system_t), abstract :: classical_particles_t
73 private
74 type(c_ptr), public :: output_handle(2)
75 type(space_t), public :: space
76 integer, public :: np
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
85 contains
86 procedure :: do_algorithmic_operation => classical_particles_do_algorithmic_operation
87 procedure :: is_tolerance_reached => classical_particles_is_tolerance_reached
88 procedure :: copy_quantities_to_interaction => classical_particles_copy_quantities_to_interaction
89 procedure :: update_interactions_start => classical_particles_update_interactions_start
90 procedure :: update_interactions_finish => classical_particles_update_interactions_finish
91 procedure :: output_start => classical_particles_output_start
92 procedure :: output_write => classical_particles_output_write
93 procedure :: output_finish => classical_particles_output_finish
94 procedure :: restart_write_data => classical_particles_restart_write_data
95 procedure :: restart_read_data => classical_particles_restart_read_data
96 procedure :: update_kinetic_energy => classical_particles_update_kinetic_energy
97 procedure :: center_of_mass => classical_particles_center_of_mass
98 procedure :: center_of_mass_vel => classical_particles_center_of_mass_vel
99 procedure :: center => classical_particles_center
100 procedure :: axis_large => classical_particles_axis_large
101 procedure :: axis_inertia => classical_particles_axis_inertia
102 end type classical_particles_t
103
104
105contains
106
107 ! ---------------------------------------------------------
112 ! ---------------------------------------------------------
113 subroutine classical_particles_init(this, np)
114 class(classical_particles_t), intent(inout) :: this
115 integer, intent(in) :: np
116
119 this%np = np
120 safe_allocate(this%mass(1:np))
121 this%mass = m_zero
122 safe_allocate(this%pos(1:this%space%dim, 1:np))
123 this%pos = m_zero
124 safe_allocate(this%vel(1:this%space%dim, 1:np))
125 this%vel = m_zero
126 safe_allocate(this%tot_force(1:this%space%dim, 1:np))
127 this%tot_force = m_zero
128 safe_allocate(this%fixed(1:np))
129 this%fixed = .false. ! By default we let the particles move.
130 safe_allocate(this%lj_epsilon(1:np))
131 this%lj_epsilon = m_zero
132 safe_allocate(this%lj_sigma(1:np))
133 this%lj_sigma = m_zero
134
135 ! No interaction directly supported by the classical particles (but classes
136 ! that extend this one can add their own)
137 allocate(this%supported_interactions(0))
138 allocate(this%supported_interactions_as_partner(0))
139
140 ! Quantities updated by the algorithm
141 call this%quantities%add(quantity_t("position", updated_on_demand = .false.))
142 call this%quantities%add(quantity_t("velocity", updated_on_demand = .false.))
143
144 ! Other quantities
145 call this%quantities%add(quantity_t("mass", updated_on_demand = .false., always_available=.true.))
146
148 end subroutine classical_particles_init
149
150 ! ---------------------------------------------------------
151 subroutine classical_particles_copy(this, cp_in)
152 class(classical_particles_t), intent(out) :: this
153 class(classical_particles_t), intent(in) :: cp_in
154
156
157 this%np = cp_in%np
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)
163
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
174 ! ---------------------------------------------------------
175 subroutine classical_particles_init_interaction(this, interaction)
176 class(classical_particles_t), target, intent(inout) :: this
177 class(interaction_surrogate_t), intent(inout) :: interaction
178
181 select type (interaction)
182 class default
183 message(1) = "Trying to initialize an unsupported interaction by classical particles."
184 call messages_fatal(1, namespace=this%namespace)
185 end select
190 ! ---------------------------------------------------------
191 logical function classical_particles_do_algorithmic_operation(this, operation, updated_quantities) result(done)
192 class(classical_particles_t), intent(inout) :: this
193 class(algorithmic_operation_t), intent(in) :: operation
194 character(len=:), allocatable, intent(out) :: updated_quantities(:)
195
196 integer :: ii, ip
197 real(real64), allocatable :: tmp_pos(:,:,:), tmp_vel(:,:,:)
198 real(real64) :: factor
199
201 call profiling_in(trim(this%namespace%get())//":"//trim(operation%id))
202
203 done = .true.
204 select type (prop => this%algo)
205 class is (propagator_t)
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)
211
212 case (verlet_start)
213 if (.not. this%prop_data%initialized) then
214 call this%prop_data%initialize(prop, this%space%dim, this%np)
215 do ip = 1, this%np
216 if (this%fixed(ip)) then
217 this%prop_data%acc(:, ip) = m_zero
218 else
219 this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
220 end if
221 end do
222 end if
223
224 case (verlet_finish)
225 call this%prop_data%end()
226
227 case (beeman_finish)
228 call this%prop_data%end()
229
230 case (verlet_update_pos)
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"]
234
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)
238 end do
239 do ip = 1, this%np
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
243 else
244 this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
245 end if
246 end do
247
248 case (verlet_compute_vel)
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"]
252
253 case (beeman_start)
254 if (.not. this%prop_data%initialized) then
255 call this%prop_data%initialize(prop, this%space%dim, this%np)
256 do ip = 1, this%np
257 if (this%fixed(ip)) then
258 this%prop_data%acc(:, ip) = m_zero
259 else
260 this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
261 end if
262 this%prop_data%prev_acc(:, ip, 1) = this%prop_data%acc(:, ip)
263 end do
264 end if
265
266 case (beeman_predict_pos)
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))
269
270 if (.not. prop%predictor_corrector) then
271 updated_quantities = ["position"]
272 end if
273
274 case (beeman_predict_vel)
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"]
279
280 case (beeman_correct_pos)
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"]
286
287 case (beeman_correct_vel)
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"]
291
292 case (expmid_2step_start)
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)
297 end if
298
300 call this%prop_data%end()
301
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"]
310
311 case (update_hamiltonian)
312 do ip = 1, this%np
313 if (this%fixed(ip)) then
314 this%prop_data%hamiltonian_elements(:, ip) = m_zero
315 else
316 this%prop_data%hamiltonian_elements(:, ip) = this%tot_force(:, ip) / (this%mass(ip) * this%pos(:, ip) + m_epsilon)
317 end if
318 end do
319
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))
323 ! apply exponential - at some point this could use the machinery of
324 ! exponential_apply (but this would require a lot of boilerplate code
325 ! like a Hamiltonian class etc)
326 ! prop_data%save_pos/vel contain the state at t - this is the state we want to
327 ! apply the Hamiltonian to
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)
332 ! compute exponential with Taylor expansion
333 factor = m_one
334 do ii = 1, 4
335 factor = factor * prop%dt / ii
336 do ip = 1, this%np
337 ! apply hamiltonian
338 tmp_pos(:, ip, 2) = tmp_vel(:, ip, 1)
339 tmp_vel(:, ip, 2) = this%prop_data%hamiltonian_elements(:, ip) * tmp_pos(:, ip, 1)
340 ! swap temporary variables
341 tmp_pos(:, ip, 1) = tmp_pos(:, ip, 2)
342 tmp_vel(:, ip, 1) = tmp_vel(:, ip, 2)
343 ! accumulate components of Taylor expansion
344 this%pos(:, ip) = this%pos(:, ip) + factor * tmp_pos(:, ip, 1)
345 this%vel(:, ip) = this%vel(:, ip) + factor * tmp_vel(:, ip, 1)
346 end do
347 end do
348 safe_deallocate_a(tmp_pos)
349 safe_deallocate_a(tmp_vel)
350 updated_quantities = ["position", "velocity"]
351
353 ! only correct for dt/2 if not converged yet
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"]
358 end if
359
360 case default
361 done = .false.
362 end select
363
364 class default
365 done = .false.
366 end select
367
368 call profiling_out(trim(this%namespace%get())//":"//trim(operation%id))
371
372 ! ---------------------------------------------------------
373 logical function classical_particles_is_tolerance_reached(this, tol) result(converged)
374 class(classical_particles_t), intent(in) :: this
375 real(real64), intent(in) :: tol
376
377 integer :: ip
378 real(real64) :: change, max_change
379
381
382 ! Here we put the criterion that maximum acceleration change is below the tolerance
383 max_change = m_zero
384 do ip = 1, this%np
385 change = sum((this%prop_data%prev_tot_force(1:this%space%dim, ip) - this%tot_force(1:this%space%dim, ip))**2) / &
386 this%mass(ip)
387 if (change > max_change) then
388 max_change = change
389 end if
390 end do
391 converged = max_change < tol**2
392
393 write(message(1), '(a, e13.6, a, e13.6)') "Debug: -- Maximum change in acceleration ", &
394 sqrt(max_change), " and tolerance ", tol
395 call messages_info(1, namespace=this%namespace, debug_only=.true.)
396
399
400 ! ---------------------------------------------------------
401 subroutine classical_particles_update_quantity(this, label)
402 class(classical_particles_t), intent(inout) :: this
403 character(len=*), intent(in) :: label
404
406
407 select case (label)
408 case default
409 message(1) = "Incompatible quantity."
410 call messages_fatal(1, namespace=this%namespace)
411 end select
412
415
416 ! ---------------------------------------------------------
417 subroutine classical_particles_init_interaction_as_partner(partner, interaction)
418 class(classical_particles_t), intent(in) :: partner
419 class(interaction_surrogate_t), intent(inout) :: interaction
420
422
423 select type (interaction)
424 class default
425 message(1) = "Unsupported interaction."
426 call messages_fatal(1, namespace=partner%namespace)
427 end select
428
431
432 ! ---------------------------------------------------------
433 subroutine classical_particles_copy_quantities_to_interaction(partner, interaction)
434 class(classical_particles_t), intent(inout) :: partner
435 class(interaction_surrogate_t), intent(inout) :: interaction
436
438
439 select type (interaction)
440 class default
441 message(1) = "Unsupported interaction."
442 call messages_fatal(1, namespace=partner%namespace)
443 end select
444
447
448 ! ---------------------------------------------------------
450 class(classical_particles_t), intent(inout) :: this
451
453
454 ! Store previous force, as it is used as SCF criterium
455 select type (prop => this%algo)
456 class is (propagator_t)
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)
459 end if
460 end select
461
464
465 ! ---------------------------------------------------------
467 class(classical_particles_t), intent(inout) :: this
468
469 type(interaction_iterator_t) :: iter
470
472
473 ! Compute the total force acting on the classical particles
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())
478 class is (force_interaction_t)
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)
481 end select
482 end do
483
486
487 ! ---------------------------------------------------------
488 subroutine classical_particles_output_start(this)
489 class(classical_particles_t), intent(inout) :: this
490
491 integer :: iteration
492 real(real64) :: dt
493
495
496 select type (algo => this%algo)
497 class is (propagator_t)
498 dt = algo%dt
499 end select
500
501 iteration = this%iteration%counter()
502 if (iteration > 0) then
503 ! add one if restarting because the output routine is only called at the end of the timestep
504 iteration = iteration + 1
505 end if
506 ! Create output handle
507 call io_mkdir('td.general', this%namespace)
508 if (mpi_grp_is_root(this%grp)) then
509 call write_iter_init(this%output_handle(output_coordinates), iteration, dt, &
510 trim(io_workpath("td.general/coordinates", this%namespace)))
511 call write_iter_init(this%output_handle(output_energy), iteration, dt, &
512 trim(io_workpath("td.general/energy", this%namespace)))
513 end if
514
515 ! Output info for first iteration
516 if (iteration == 0) then
517 call this%output_write()
518 end if
519
522
523 ! ---------------------------------------------------------
524 subroutine classical_particles_output_finish(this)
525 class(classical_particles_t), intent(inout) :: this
528
529 if (mpi_grp_is_root(this%grp)) then
530 call write_iter_end(this%output_handle(output_coordinates))
531 call write_iter_end(this%output_handle(output_energy))
532 end if
533
536
537 ! ---------------------------------------------------------
538 subroutine classical_particles_output_write(this)
539 class(classical_particles_t), intent(inout) :: this
540
541 integer :: idir, ii, ip
542 character(len=50) :: aux
543 real(real64) :: tmp(this%space%dim)
544
545 if (.not. mpi_grp_is_root(this%grp)) return ! only first node outputs
546
548
549 if (this%iteration%counter() == 0) then
550 ! header
551 call write_iter_clear(this%output_handle(output_coordinates))
552 call write_iter_string(this%output_handle(output_coordinates), &
553 '############################################################################')
554 call write_iter_nl(this%output_handle(output_coordinates))
555 call write_iter_string(this%output_handle(output_coordinates),'# HEADER')
556 call write_iter_nl(this%output_handle(output_coordinates))
557
558 ! first line: column names
559 call write_iter_header_start(this%output_handle(output_coordinates))
560 do ip = 1, this%np
561 do idir = 1, this%space%dim
562 write(aux, '(a2,i3,a1,i3,a1)') 'x(', ip, ',', idir, ')'
563 call write_iter_header(this%output_handle(output_coordinates), aux)
564 end do
565 end do
566 do ip = 1, this%np
567 do idir = 1, this%space%dim
568 write(aux, '(a2,i3,a1,i3,a1)') 'v(', ip, ',', idir, ')'
569 call write_iter_header(this%output_handle(output_coordinates), aux)
570 end do
571 end do
572 do ip = 1, this%np
573 do idir = 1, this%space%dim
574 write(aux, '(a2,i3,a1,i3,a1)') 'f(', ip, ',', idir, ')'
575 call write_iter_header(this%output_handle(output_coordinates), aux)
576 end do
577 end do
578 call write_iter_nl(this%output_handle(output_coordinates))
579
580 ! second line: units
581 call write_iter_string(this%output_handle(output_coordinates), '#[Iter n.]')
582 call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%time)) // ']')
583 do ip = 1, this%np
584 do idir = 1, this%space%dim
585 call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%length)) // ']')
586 end do
587 end do
588 do ip = 1, this%np
589 do idir = 1, this%space%dim
590 call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%velocity)) // ']')
591 end do
592 end do
593 do ip = 1, this%np
594 do idir = 1, this%space%dim
595 call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%force)) // ']')
596 end do
597 end do
598 call write_iter_nl(this%output_handle(output_coordinates))
599
600 call write_iter_string(this%output_handle(output_coordinates), &
601 '############################################################################')
602 call write_iter_nl(this%output_handle(output_coordinates))
603 end if
604
605 call write_iter_start(this%output_handle(output_coordinates))
606
607 ! Position
608 do ip = 1, this%np
609 tmp(1:this%space%dim) = units_from_atomic(units_out%length, this%pos(1:this%space%dim, ip))
610 call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
611 end do
612 ! Velocity
613 do ip = 1, this%np
614 tmp(1:this%space%dim) = units_from_atomic(units_out%velocity, this%vel(1:this%space%dim, ip))
615 call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
616 end do
617 ! Force
618 do ip = 1, this%np
619 tmp(1:this%space%dim) = units_from_atomic(units_out%force, this%tot_force(1:this%space%dim, ip))
620 call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
621 end do
622
623 call write_iter_nl(this%output_handle(output_coordinates))
624
625 ! Energies
626 if (this%iteration%counter() == 0) then
627 ! header
628 call write_iter_clear(this%output_handle(output_energy))
629 call write_iter_string(this%output_handle(output_energy), &
630 '############################################################################')
631 call write_iter_nl(this%output_handle(output_energy))
632 call write_iter_string(this%output_handle(output_energy),'# HEADER')
633 call write_iter_nl(this%output_handle(output_energy))
634
635 ! first line: column names
636 call write_iter_header_start(this%output_handle(output_energy))
637 call write_iter_header(this%output_handle(output_energy), 'Total')
638 call write_iter_header(this%output_handle(output_energy), 'Kinetic')
639 call write_iter_header(this%output_handle(output_energy), 'Potential')
640 call write_iter_header(this%output_handle(output_energy), 'Internal')
641 call write_iter_nl(this%output_handle(output_energy))
642
643 ! second line: units
644 call write_iter_string(this%output_handle(output_energy), '#[Iter n.]')
645 call write_iter_header(this%output_handle(output_energy), '[' // trim(units_abbrev(units_out%time)) // ']')
646 do ii = 1, 4
647 call write_iter_header(this%output_handle(output_energy), '[' // trim(units_abbrev(units_out%energy)) // ']')
648 end do
649 call write_iter_nl(this%output_handle(output_energy))
650
651 call write_iter_string(this%output_handle(output_energy), &
652 '############################################################################')
653 call write_iter_nl(this%output_handle(output_energy))
654 end if
655 call write_iter_start(this%output_handle(output_energy))
656
657 call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%total_energy), 1)
658 call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%kinetic_energy), 1)
659 call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%potential_energy), 1)
660 call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%internal_energy), 1)
661 call write_iter_nl(this%output_handle(output_energy))
662
665
666 ! ---------------------------------------------------------
668 class(classical_particles_t), intent(inout) :: this
669
670 integer :: restart_file_unit
671
673
674 if (mpi_grp_is_root(this%grp)) then
675 call write_iter_flush(this%output_handle(output_coordinates))
676 call write_iter_flush(this%output_handle(output_energy))
677 end if
678
679 call profiling_in(trim(this%namespace%get())//":"//"RESTART_WRITE")
680
681 call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
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(:,:)
688 call io_close(restart_file_unit)
689
690 if (this%iteration%counter() > 0) then
691 ! only initialized after the first time step
692 call this%prop_data%restart_write(this%namespace, this%algo)
693 end if
694
695 message(1) = "Successfully wrote restart data for system "//trim(this%namespace%get())
696 call messages_info(1, namespace=this%namespace)
697
698 call profiling_out(trim(this%namespace%get())//":"//"RESTART_WRITE")
701
702 ! ---------------------------------------------------------
703 logical function classical_particles_restart_read_data(this)
704 class(classical_particles_t), intent(inout) :: this
705
706 integer :: restart_file_unit
707
709 call profiling_in(trim(this%namespace%get())//":"//"RESTART_READ")
710
711 call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
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(:,:)
719 call io_close(restart_file_unit)
720 call this%prop_data%initialize(this%algo, this%space%dim, this%np)
721 if (this%iteration%counter() > 0) then
722 ! only initialized after the first time step
723 classical_particles_restart_read_data = this%prop_data%restart_read(this%namespace, this%algo)
724 else
726 end if
727 else
728 ! could not open file
730 end if
731
733 message(1) = "Successfully read restart data for system "//trim(this%namespace%get())
734 ! namespace should be added here at some point
735 call messages_info(1)
736 else
737 message(1) = "Failed to read restart data for system "//trim(this%namespace%get())
738 ! namespace should be added here at some point
739 call messages_info(1)
740 end if
741
742 call profiling_out(trim(this%namespace%get())//":"//"RESTART_READ")
745
746 ! ---------------------------------------------------------
748 class(classical_particles_t), intent(inout) :: this
749
750 integer :: ip
751
753
754 this%kinetic_energy = m_zero
755 do ip = 1, this%np
756 this%kinetic_energy = this%kinetic_energy + m_half * this%mass(ip) * sum(this%vel(:, ip)**2)
757 end do
758
761
762 ! ---------------------------------------------------------
763 function classical_particles_center_of_mass(this, mask, pseudo) result(pos)
764 class(classical_particles_t), intent(in) :: this
765 logical, optional, intent(in) :: mask(:)
766 logical, optional, intent(in) :: pseudo
767 real(real64) :: pos(this%space%dim)
768
769 real(real64) :: mass, total_mass
770 integer :: ip
771
773
774 assert(.not. this%space%is_periodic())
775
776 pos = m_zero
777 total_mass = m_zero
778 mass = m_one
779 do ip = 1, this%np
780 if (present(mask)) then
781 if (.not. mask(ip)) cycle
782 end if
783 if (.not. optional_default(pseudo, .false.)) then
784 mass = this%mass(ip)
785 end if
786 pos = pos + mass*this%pos(:, ip)
787 total_mass = total_mass + mass
788 end do
789 pos = pos/total_mass
790
793
794 ! ---------------------------------------------------------
795 function classical_particles_center_of_mass_vel(this) result(vel)
796 class(classical_particles_t), intent(in) :: this
797 real(real64) :: vel(this%space%dim)
798
799 real(real64) :: mass, total_mass
800 integer :: ip
801
803
804 vel = m_zero
805 total_mass = m_zero
806 do ip = 1, this%np
807 mass = this%mass(ip)
808 total_mass = total_mass + mass
809 vel = vel + mass*this%vel(:, ip)
810 end do
811 vel = vel/total_mass
812
815
816 ! ---------------------------------------------------------
817 function classical_particles_center(this) result(pos)
818 class(classical_particles_t), intent(in) :: this
819 real(real64) :: pos(this%space%dim)
820
821 real(real64) :: xmin(this%space%dim), xmax(this%space%dim)
822 integer :: ip, idir
823
825
826 xmin = m_huge
827 xmax = -m_huge
828 do ip = 1, this%np
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)
832 end do
833 end do
834
835 pos = (xmax + xmin)/m_two
836
838 end function classical_particles_center
839
840 ! ---------------------------------------------------------
841 subroutine classical_particles_axis_large(this, x, x2)
842 class(classical_particles_t), intent(in) :: this
843 real(real64), intent(out) :: x(this%space%dim)
844 real(real64), intent(out) :: x2(this%space%dim)
845
846 integer :: ip, jp
847 real(real64) :: rmax, r, r2
848
850
851 assert(.not. this%space%is_periodic())
852
853 ! first get the further apart atoms
854 rmax = -m_huge
855 do ip = 1, this%np
856 do jp = 1, this%np/2 + 1
857 r = norm2(this%pos(:, ip) - this%pos(:, jp))
858 if (r > rmax) then
859 rmax = r
860 x = this%pos(:, ip) - this%pos(:, jp)
861 end if
862 end do
863 end do
864 x = x /norm2(x)
865
866 ! now let us find out what is the second most important axis
867 rmax = -m_huge
868 do ip = 1, this%np
869 r2 = sum(x * this%pos(:, ip))
870 r = norm2(this%pos(:, ip) - r2*x)
871 if (r > rmax) then
872 rmax = r
873 x2 = this%pos(:, ip) - r2*x
874 end if
875 end do
876
878 end subroutine classical_particles_axis_large
879
880 ! ---------------------------------------------------------
883 subroutine classical_particles_axis_inertia(this, x, x2, pseudo)
884 class(classical_particles_t), intent(in) :: this
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
891 type(unit_t) :: unit
892
894
895 assert(.not. this%space%is_periodic())
896
897 ! first calculate the inertia tensor
898 tinertia = m_zero
899 mass = m_one
900 do ip = 1, this%np
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)
905 end do
906 end do
907
908 unit = units_out%length**2
909 ! note: we always use amu for atomic masses, so no unit conversion to/from atomic is needed.
910 if (pseudo) then
911 write(message(1),'(a)') 'Moment of pseudo-inertia tensor [' // trim(units_abbrev(unit)) // ']'
912 else
913 write(message(1),'(a)') 'Moment of inertia tensor [amu*' // trim(units_abbrev(unit)) // ']'
914 end if
915 call messages_info(1, namespace=this%namespace)
916 call output_tensor(tinertia, this%space%dim, unit, write_average = .true., namespace=this%namespace)
917
918 call lalg_eigensolve(this%space%dim, tinertia, eigenvalues)
919
920 write(message(1),'(a,6f25.6)') 'Eigenvalues: ', units_from_atomic(unit, eigenvalues)
921 call messages_info(1, namespace=this%namespace)
922
923 ! make a choice to fix the sign of the axis.
924 do ii = 1, 2
925 jj = maxloc(abs(tinertia(:,ii)), dim = 1)
926 if (tinertia(jj,ii) < m_zero) tinertia(:,ii) = -tinertia(:,ii)
927 end do
928 x = tinertia(:,1)
929 x2 = tinertia(:,2)
930
933
934 ! ---------------------------------------------------------
935 subroutine classical_particles_end(this)
936 class(classical_particles_t), intent(inout) :: this
937
939
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)
947
948 call system_end(this)
949
951 end subroutine classical_particles_end
952
954
955!! Local Variables:
956!! mode: f90
957!! coding: utf-8
958!! End:
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
Definition: write_iter.F90:167
This module implements the basic elements defining algorithms.
Definition: algorithm.F90:141
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)
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
Definition: global.F90:190
real(real64), parameter, public m_huge
Definition: global.F90:206
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_epsilon
Definition: global.F90:204
character(len= *), parameter, public td_dir
Definition: global.F90:251
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
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
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
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
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_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
This module implements the multisystem debug functionality.
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
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.
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
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
subroutine, public output_tensor(tensor, ndim, unit, write_average, iunit, namespace)
Definition: utils.F90:242
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
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.
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)