-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.F90
126 lines (92 loc) · 5.08 KB
/
main.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
program main
use, intrinsic :: iso_fortran_env, only: DP => REAL64, I64 => INT64
use :: therm_m, only: calc_KE, calc_vdw_pbc
use :: init_cond_m, only: init_positions_sc, bimodal_dist_velocities
use :: integrtors_m, only: velocity_verlet, euler_integrator
use :: writers_m, only: write_velocities
use :: interface_m, only: databloc_params_t
use :: readers_m, only: read_nml
implicit none
! ~ Memory definition ~
! Derived types
type(databloc_params_t) :: datablock
! Scalar values
real(kind=DP) :: time0, time1
integer(kind=I64) :: log_unit, vel_unit, rdf_unit, msd_unit, frame_unit, gen_log_unit
integer :: status_cli
character(len=2048) :: log_name, vel_name, rdf_name, msd_name, frame_name, nml_path, gen_log_name
! Arrays
real(kind=DP), allocatable, dimension(:, :) :: positions, velocities, init_position
call get_command_argument(1, nml_path, status=status_cli)
if (status_cli /= 0) stop 1
! Inicialitzem el RNG
call random_seed()
! ~ Definim parametres i inicialitzem parametres ~
call read_nml(param_file=nml_path, datablock=datablock)
gen_log_name = trim(datablock%sim_name) // "_general.log"
! ----------------------------------------------------------------------------------------------------------------------------------------
! ~ Realitzem l'equilibrat del sistema ~
datablock%n_steps = 100000_I64
datablock%ref_temp = 100.0_dp
log_name = trim(datablock%sim_name) // "_heating.log"
vel_name = trim(datablock%sim_name) // "_heating.vel"
! Files
open(newunit=gen_log_unit, file=trim(gen_log_name), access='sequential', action='write',&
status='replace', form='formatted')
open(newunit=log_unit, file=trim(log_name), access='sequential', action='write',&
status='replace', form='formatted')
open(newunit=vel_unit, file=trim(vel_name), access='sequential', action='write',&
status='replace', form='formatted')
! ~ Inicialitzem les posicions i velocitats~
allocate(positions(3, datablock%n_particles), velocities(3, datablock%n_particles))
allocate(init_position(3, datablock%n_particles))
call init_positions_sc(rho=datablock%density, pos=positions)
call bimodal_dist_velocities(vel=velocities, temp=datablock%ref_temp)
Write(unit=gen_log_unit, fmt='(A,F16.8,A,F16.8)') "Initial potential energy=", &
calc_vdw_pbc(pos=positions, cutoff=datablock%cutoff_set, boundary=datablock%box), &
" Initial kinetic energy=", calc_KE(velocities)
call write_velocities(vel=velocities, unit_nr=vel_unit, step=0_i64)
call cpu_time(time0)
call velocity_verlet(vel=velocities, pos=positions, parambox=datablock, &
log_unit=log_unit, thermost=.TRUE.)
call cpu_time(time1)
call write_velocities(velocities, vel_unit, datablock%n_steps)
write(unit=gen_log_unit, fmt='(A,F12.8,A,F12.8)') "Execution time for heating: ", time1 - time0, " time/iteration: ", (time1 - time0) / datablock%n_steps
close(log_unit)
close(vel_unit)
! ----------------------------------------------------------------------------------------------------------------------------------------
! ~ Realitzem la producció del sistema ~
datablock%ref_temp = 1.2_dp
datablock%n_steps = 500000_i64
init_position = positions
! Files
log_name = trim(datablock%sim_name) // "_prod.log"
vel_name = trim(datablock%sim_name) // "_prod.vel"
frame_name = trim(datablock%sim_name) // "_prod.xyz"
rdf_name = trim(datablock%sim_name) // "_rdf.dat"
msd_name = trim(datablock%sim_name) // "_msd.dat"
open(newunit=log_unit, file=trim(log_name), access='sequential', action='write',&
status='replace', form='formatted')
open(newunit=vel_unit, file=trim(vel_name), access='sequential', action='write',&
status='replace', form='formatted')
open(newunit=rdf_unit, file=trim(rdf_name), access='sequential', action='write',&
status='replace', form='formatted')
open(newunit=msd_unit, file=trim(msd_name), access='sequential', action='write',&
status='replace', form='formatted')
open(newunit=frame_unit, file=trim(frame_name), access='sequential', action='write',&
status='replace', form='formatted')
call cpu_time(time0)
call velocity_verlet(vel=velocities, pos=positions, parambox=datablock, &
log_unit=log_unit, thermost=.TRUE., rdf_unit=rdf_unit, &
msd_unit=msd_unit, init_pos=init_position, frame_unit=frame_unit)
call cpu_time(time1)
write(unit=gen_log_unit, fmt='(A,F12.8,A,F12.8)') "Execution time for production: ", time1 - time0, " time/iteration: ", (time1 - time0) / datablock%n_steps
close(log_unit)
close(vel_unit)
close(rdf_unit)
close(msd_unit)
close(frame_unit)
close(gen_log_unit)
deallocate(positions)
deallocate(velocities)
end program main