1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief module handles definition of the tree nodes for the global and
8!>      the subtrees binary tree
9!>                   parent element
10!>                      /      \
11!>      accepted (acc) /        \  not accepted (nacc)
12!>                    /          \
13!>                  child       child
14!>                   / \         / \
15!>
16!>      tree creation assuming acceptance (acc) AND rejectance (nacc)
17!>        of configuration
18!>      if configuration is accepted: new configuration (child on acc) on basis
19!>        of last configuration (one level up)
20!>      if configuration is rejected: child on nacc on basis of last accepted
21!>        element (last element which is on acc brach of its parent element)
22!>      The global tree handles all configurations of different subtrees.
23!>      The structure element "conf" is an array related to the temperature
24!>        (sorted) and points to the subtree elements.
25!> \par History
26!>      11.2012 created [Mandes Schoenherr]
27!> \author Mandes
28! **************************************************************************************************
29
30MODULE tmc_tree_types
31   USE kinds,                           ONLY: dp
32#include "../base/base_uses.f90"
33
34   IMPLICIT NONE
35
36   PRIVATE
37
38   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_tree_types'
39
40   PUBLIC :: tree_type, global_tree_type
41   PUBLIC :: elem_list_type, elem_array_type, gt_elem_list_type
42   PUBLIC :: add_to_list, clean_list
43   PUBLIC :: read_subtree_elem_unformated, write_subtree_elem_unformated
44
45   !-- tree element status
46   INTEGER, PARAMETER, PUBLIC :: status_created = 100
47   INTEGER, PARAMETER, PUBLIC :: status_calculate_energy = 101
48   INTEGER, PARAMETER, PUBLIC :: status_calc_approx_ener = 102
49
50   INTEGER, PARAMETER, PUBLIC :: status_calculate_NMC_steps = 111
51   INTEGER, PARAMETER, PUBLIC :: status_calculate_MD = 112
52   INTEGER, PARAMETER, PUBLIC :: status_calculated = 113
53
54   INTEGER, PARAMETER, PUBLIC :: status_accepted_result = 123
55   INTEGER, PARAMETER, PUBLIC :: status_accepted = 122
56   INTEGER, PARAMETER, PUBLIC :: status_rejected = 121
57   INTEGER, PARAMETER, PUBLIC :: status_rejected_result = 120
58
59   INTEGER, PARAMETER, PUBLIC :: status_cancel_nmc = 133
60   INTEGER, PARAMETER, PUBLIC :: status_cancel_ener = 132
61   INTEGER, PARAMETER, PUBLIC :: status_canceled_nmc = 131
62   INTEGER, PARAMETER, PUBLIC :: status_canceled_ener = 130
63
64   INTEGER, PARAMETER, PUBLIC :: status_deleted = 140
65   INTEGER, PARAMETER, PUBLIC :: status_deleted_result = 141
66
67   !-- dimension status (for e.g. dividing atoms in sub box)
68   INTEGER, PARAMETER, PUBLIC :: status_ok = 42
69   INTEGER, PARAMETER, PUBLIC :: status_frozen = -1
70   INTEGER, PARAMETER, PUBLIC :: status_proton_disorder = 1
71
72   !-- subtree element
73   TYPE tree_type
74      TYPE(tree_type), POINTER                :: parent => NULL() ! points to element one level up
75      !-- acc..accepted goes to next level (next step),
76      !   nacc..not accepted takes an alternative configutation
77      TYPE(tree_type), POINTER                :: acc => NULL(), nacc => NULL()
78      !-- type of MC move (swap is handled only in global tree)
79      INTEGER                                  :: move_type
80      !-- status (e.g. calculated, MD calculation, accepted...)
81      INTEGER                                  :: stat = status_created
82      REAL(KIND=dp), DIMENSION(:), POINTER     :: subbox_center
83      REAL(KIND=dp), DIMENSION(:), POINTER     :: pos ! position array
84      INTEGER, DIMENSION(:), POINTER           :: mol ! specifies the molecules the atoms participate
85      REAL(KIND=dp), DIMENSION(:), POINTER     :: vel ! velocity array
86      REAL(KIND=dp), DIMENSION(:), POINTER     :: frc ! force array
87      REAL(KIND=dp), DIMENSION(:), POINTER     :: dipole ! dipole moments array
88      INTEGER, DIMENSION(:), POINTER           :: elem_stat ! status for every dimension
89      INTEGER                                  :: nr ! tree node number
90      REAL(KIND=dp), DIMENSION(3, 2, 3)        :: rng_seed ! random seed for childs
91      !-- remembers which subtree number element is from
92      INTEGER                                  :: sub_tree_nr
93      !-- remembers the temperature the configurational change (NMC) is done with
94      INTEGER                                  :: temp_created = 0
95      !-- pointer to counter of next subtree element number
96      INTEGER, POINTER                         :: next_elem_nr
97      !-- for calculating the NPT ensamble, variable box sizes are necessary.
98      REAL(KIND=dp), DIMENSION(:), POINTER     :: box_scale
99      REAL(KIND=dp)                            :: potential ! potential energy
100      !-- potential energy calculated using (MD potential) cp2k input file
101      REAL(KIND=dp)                            :: e_pot_approx = 0.0_dp
102      !-- kinetic energy (espacially for HMC, where the velocities are respected)
103      REAL(KIND=dp)                            :: ekin
104      !-- kinetic energy before md steps (after gaussian velocity change)
105      REAL(KIND=dp)                            :: ekin_before_md
106      !-- estimated energies are stored in loop order in this array
107      REAL(KIND=dp), DIMENSION(4)              :: scf_energies
108      !-- counter to get last position in the array loop
109      INTEGER                                  :: scf_energies_count
110      !-- list of global tree elements referint to that node (reference back to global tree)
111      !   if no reference exist anymore, global tree element can be deleted
112      TYPE(gt_elem_list_type), POINTER         :: gt_nodes_references => NULL()
113   END TYPE tree_type
114
115   ! type for global tree element list in tree elements
116   TYPE gt_elem_list_type
117      TYPE(global_tree_type), POINTER         :: gt_elem => NULL()
118      TYPE(gt_elem_list_type), POINTER        :: next => NULL()
119   END TYPE gt_elem_list_type
120
121   TYPE elem_list_type
122      TYPE(tree_type), POINTER      :: elem => NULL()
123      TYPE(elem_list_type), POINTER :: next => NULL()
124      INTEGER                        :: temp_ind
125      INTEGER                        :: nr
126   END TYPE elem_list_type
127
128   !-- array with subtree elements
129   TYPE elem_array_type
130      TYPE(tree_type), POINTER :: elem => NULL()
131      LOGICAL                   :: busy = .FALSE.
132      LOGICAL                   :: canceled = .FALSE.
133      REAL(KIND=dp)             :: start_time
134   END TYPE elem_array_type
135
136   !-- global tree element
137   TYPE global_tree_type
138      TYPE(global_tree_type), POINTER :: parent => NULL() ! points to element one level up
139      !-- acc..accepted goes to next level (next step),
140      !   nacc..not accepted takes an alternative configutation
141      TYPE(global_tree_type), POINTER :: acc => NULL(), nacc => NULL()
142      !-- status (e.g. calculated, MD calculation, accepted...)
143      INTEGER                                      :: stat = -99
144      !-- remember if configuration in node are swaped
145      LOGICAL                                      :: swaped
146      !-- stores the index of the configuration (temperature)
147      !   which is changed
148      INTEGER                                      :: mv_conf = -54321
149      !-- stores the index of the configuration (temp.) which should change next
150      INTEGER                                      :: mv_next_conf = -2345
151      !-- list of pointes to subtree elements (Temp sorting)
152      TYPE(elem_array_type), DIMENSION(:), ALLOCATABLE :: conf
153      !-- remembers if last configuration is assumed to be accepted or rejected (next branc in tree);
154      !   In case of swaping, it shows if the configuration of a certain temperature is assumed
155      !   to be acc/rej (which branch is followed at the last modification of the conf of this temp.
156      !TODO store conf_n_acc in a bitshifted array to decrease the size (1Logical = 1Byte)
157      LOGICAL, DIMENSION(:), ALLOCATABLE           :: conf_n_acc
158      INTEGER :: nr ! tree node number
159      REAL(KIND=dp), DIMENSION(3, 2, 3)            :: rng_seed ! random seed for childs
160      !-- random number for acceptance check
161      REAL(KIND=dp)                                :: rnd_nr
162      !-- approximate probability of acceptance will be adapted while calculating the exact energy
163      REAL(KIND=dp)                                :: prob_acc ! estimated acceptance probability
164      REAL(KIND=dp)                                :: Temp ! temperature for simulated annealing
165   END TYPE global_tree_type
166
167CONTAINS
168
169! **************************************************************************************************
170!> \brief add a certain element to the specified element list at the beginning
171!> \param elem the sub tree element, to be added
172!> \param list  ...
173!> \param temp_ind ...
174!> \param nr ...
175!> \author Mandes 11.2012
176! **************************************************************************************************
177   SUBROUTINE add_to_list(elem, list, temp_ind, nr)
178      TYPE(tree_type), POINTER                           :: elem
179      TYPE(elem_list_type), POINTER                      :: list
180      INTEGER, OPTIONAL                                  :: temp_ind, nr
181
182      CHARACTER(LEN=*), PARAMETER :: routineN = 'add_to_list', routineP = moduleN//':'//routineN
183
184      TYPE(elem_list_type), POINTER                      :: last, list_elem_tmp
185
186      NULLIFY (list_elem_tmp, last)
187
188      CPASSERT(ASSOCIATED(elem))
189
190      ALLOCATE (list_elem_tmp)
191      list_elem_tmp%elem => elem
192      list_elem_tmp%next => NULL()
193      IF (PRESENT(temp_ind)) THEN
194         list_elem_tmp%temp_ind = temp_ind
195      ELSE
196         list_elem_tmp%temp_ind = -1
197      END IF
198
199      IF (PRESENT(nr)) THEN
200         list_elem_tmp%nr = nr
201      ELSE
202         list_elem_tmp%nr = -1
203      END IF
204
205      IF (ASSOCIATED(list) .EQV. .FALSE.) THEN
206         list => list_elem_tmp
207      ELSE
208         last => list
209         DO WHILE (ASSOCIATED(last%next))
210            last => last%next
211         END DO
212         last%next => list_elem_tmp
213      END IF
214
215   END SUBROUTINE add_to_list
216
217! **************************************************************************************************
218!> \brief clean a certain element element list
219!> \param list  ...
220!> \author Mandes 11.2012
221! **************************************************************************************************
222   SUBROUTINE clean_list(list)
223      TYPE(elem_list_type), POINTER                      :: list
224
225      CHARACTER(LEN=*), PARAMETER :: routineN = 'clean_list', routineP = moduleN//':'//routineN
226
227      TYPE(elem_list_type), POINTER                      :: list_elem_tmp
228
229      NULLIFY (list_elem_tmp)
230
231      DO WHILE (ASSOCIATED(list))
232         list_elem_tmp => list%next
233         DEALLOCATE (list)
234         list => list_elem_tmp
235      END DO
236   END SUBROUTINE clean_list
237
238! **************************************************************************************************
239!> \brief prints out the TMC sub tree structure element unformated in file
240!> \param elem ...
241!> \param io_unit ...
242!> \param
243!> \author Mandes 11.2012
244! **************************************************************************************************
245   SUBROUTINE write_subtree_elem_unformated(elem, io_unit)
246      TYPE(tree_type), POINTER                           :: elem
247      INTEGER                                            :: io_unit
248
249      CHARACTER(LEN=*), PARAMETER :: routineN = 'write_subtree_elem_unformated', &
250         routineP = moduleN//':'//routineN
251
252      CPASSERT(ASSOCIATED(elem))
253      CPASSERT(io_unit .GT. 0)
254      WRITE (io_unit) elem%nr, &
255         elem%sub_tree_nr, &
256         elem%stat, &
257         elem%rng_seed, &
258         elem%move_type, &
259         elem%temp_created, &
260         elem%potential, &
261         elem%e_pot_approx, &
262         elem%ekin, &
263         elem%ekin_before_md
264      CALL write_subtree_elem_darray(elem%pos, io_unit)
265      CALL write_subtree_elem_darray(elem%vel, io_unit)
266      CALL write_subtree_elem_darray(elem%frc, io_unit)
267      CALL write_subtree_elem_darray(elem%box_scale, io_unit)
268      CALL write_subtree_elem_darray(elem%dipole, io_unit)
269   END SUBROUTINE write_subtree_elem_unformated
270
271! **************************************************************************************************
272!> \brief reads the TMC sub tree structure element unformated in file
273!> \param elem ...
274!> \param io_unit ...
275!> \param
276!> \author Mandes 11.2012
277! **************************************************************************************************
278   SUBROUTINE read_subtree_elem_unformated(elem, io_unit)
279      TYPE(tree_type), POINTER                           :: elem
280      INTEGER                                            :: io_unit
281
282      CHARACTER(LEN=*), PARAMETER :: routineN = 'read_subtree_elem_unformated', &
283         routineP = moduleN//':'//routineN
284
285      CPASSERT(ASSOCIATED(elem))
286      CPASSERT(io_unit .GT. 0)
287
288      READ (io_unit) elem%nr, &
289         elem%sub_tree_nr, &
290         elem%stat, &
291         elem%rng_seed, &
292         elem%move_type, &
293         elem%temp_created, &
294         elem%potential, &
295         elem%e_pot_approx, &
296         elem%ekin, &
297         elem%ekin_before_md
298      CALL read_subtree_elem_darray(elem%pos, io_unit)
299      CALL read_subtree_elem_darray(elem%vel, io_unit)
300      CALL read_subtree_elem_darray(elem%frc, io_unit)
301      CALL read_subtree_elem_darray(elem%box_scale, io_unit)
302      CALL read_subtree_elem_darray(elem%dipole, io_unit)
303   END SUBROUTINE read_subtree_elem_unformated
304
305! **************************************************************************************************
306!> \brief ...
307!> \param array ...
308!> \param io_unit ...
309! **************************************************************************************************
310   SUBROUTINE write_subtree_elem_darray(array, io_unit)
311      REAL(KIND=dp), DIMENSION(:), POINTER               :: array
312      INTEGER                                            :: io_unit
313
314      CHARACTER(LEN=*), PARAMETER :: routineN = 'write_subtree_elem_darray', &
315         routineP = moduleN//':'//routineN
316
317      WRITE (io_unit) ASSOCIATED(array)
318      IF (ASSOCIATED(array)) THEN
319         WRITE (io_unit) SIZE(array)
320         WRITE (io_unit) array
321      END IF
322   END SUBROUTINE write_subtree_elem_darray
323
324! **************************************************************************************************
325!> \brief ...
326!> \param array ...
327!> \param io_unit ...
328! **************************************************************************************************
329   SUBROUTINE read_subtree_elem_darray(array, io_unit)
330      REAL(KIND=dp), DIMENSION(:), POINTER               :: array
331      INTEGER                                            :: io_unit
332
333      CHARACTER(LEN=*), PARAMETER :: routineN = 'read_subtree_elem_darray', &
334         routineP = moduleN//':'//routineN
335
336      INTEGER                                            :: i_tmp
337      LOGICAL                                            :: l_tmp
338
339      READ (io_unit) l_tmp
340      IF (l_tmp) THEN
341         READ (io_unit) i_tmp
342         IF (ASSOCIATED(array)) THEN
343            CPASSERT(SIZE(array) .EQ. i_tmp)
344         ELSE
345            ALLOCATE (array(i_tmp))
346         END IF
347         READ (io_unit) array
348      END IF
349   END SUBROUTINE read_subtree_elem_darray
350
351END MODULE tmc_tree_types
352