Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief contains miscellaneous subroutines used in the Monte Carlo runs,mostly
10 : !> geared towards changes in coordinates
11 : !> \author MJM
12 : ! **************************************************************************************************
13 : MODULE mc_coordinates
14 : USE cell_types, ONLY: cell_type,&
15 : get_cell
16 : USE cp_subsys_types, ONLY: cp_subsys_get,&
17 : cp_subsys_type
18 : USE force_env_methods, ONLY: force_env_calc_energy_force
19 : USE force_env_types, ONLY: force_env_get,&
20 : force_env_type
21 : USE kinds, ONLY: dp
22 : USE mathconstants, ONLY: pi
23 : USE mc_types, ONLY: get_mc_molecule_info,&
24 : get_mc_par,&
25 : mc_molecule_info_type,&
26 : mc_simpar_type
27 : USE message_passing, ONLY: mp_comm_type
28 : USE molecule_types, ONLY: molecule_type
29 : USE parallel_rng_types, ONLY: rng_stream_type
30 : USE particle_list_types, ONLY: particle_list_type
31 : USE physcon, ONLY: angstrom
32 : #include "../../base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : PRIVATE :: generate_avbmc_insertion
39 :
40 : PUBLIC :: generate_cbmc_swap_config, &
41 : get_center_of_mass, mc_coordinate_fold, &
42 : find_mc_test_molecule, &
43 : create_discrete_array, &
44 : check_for_overlap, &
45 : rotate_molecule, &
46 : cluster_search
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_coordinates'
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief looks for overlaps (intermolecular distances less than rmin)
54 : !> \param force_env the force environment containing the coordinates
55 : !> \param nchains the number of molecules of each type in the box
56 : !> \param nunits the number of interaction sites for each molecule
57 : !> \param loverlap returns .TRUE. if atoms overlap
58 : !> \param mol_type an array that indicates the type of each molecule
59 : !> \param cell_length the length of the box...if none is specified,
60 : !> it uses the cell found in the force_env
61 : !> \param molecule_number if present, just look for overlaps with this
62 : !> molecule
63 : !>
64 : !> Suitable for parallel use.
65 : !> \author MJM
66 : ! **************************************************************************************************
67 8650 : SUBROUTINE check_for_overlap(force_env, nchains, nunits, loverlap, mol_type, &
68 : cell_length, molecule_number)
69 :
70 : TYPE(force_env_type), POINTER :: force_env
71 : INTEGER, DIMENSION(:), INTENT(IN) :: nchains, nunits
72 : LOGICAL, INTENT(OUT) :: loverlap
73 : INTEGER, DIMENSION(:), INTENT(IN) :: mol_type
74 : REAL(KIND=dp), DIMENSION(1:3), INTENT(IN), &
75 : OPTIONAL :: cell_length
76 : INTEGER, INTENT(IN), OPTIONAL :: molecule_number
77 :
78 : CHARACTER(len=*), PARAMETER :: routineN = 'check_for_overlap'
79 :
80 : INTEGER :: handle, imol, iunit, jmol, jstart, &
81 : junit, nend, nstart, nunit, nunits_i, &
82 : nunits_j
83 : LOGICAL :: lall
84 : REAL(KIND=dp) :: dist, rmin
85 8650 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r
86 : REAL(KIND=dp), DIMENSION(1:3) :: abc, box_length, RIJ
87 : TYPE(cell_type), POINTER :: cell
88 : TYPE(cp_subsys_type), POINTER :: oldsys
89 : TYPE(particle_list_type), POINTER :: particles
90 :
91 : ! begin the timing of the subroutine
92 :
93 8650 : CALL timeset(routineN, handle)
94 :
95 8650 : NULLIFY (oldsys, particles)
96 :
97 : ! initialize some stuff
98 8650 : loverlap = .FALSE.
99 8650 : rmin = 3.57106767_dp ! 1 angstrom squared
100 :
101 : ! get the particle coordinates and the cell length
102 8650 : CALL force_env_get(force_env, cell=cell, subsys=oldsys)
103 8650 : CALL get_cell(cell, abc=abc)
104 8650 : CALL cp_subsys_get(oldsys, particles=particles)
105 :
106 57316 : ALLOCATE (r(1:3, 1:MAXVAL(nunits), 1:SUM(nchains)))
107 :
108 8650 : IF (PRESENT(cell_length)) THEN
109 82 : box_length(1:3) = cell_length(1:3)
110 : ELSE
111 8568 : box_length(1:3) = abc(1:3)
112 : END IF
113 :
114 : ! put the coordinates into an easier matrix to manipulate
115 8650 : junit = 0
116 66446 : DO imol = 1, SUM(nchains)
117 46438 : nunit = nunits(mol_type(imol))
118 162532 : DO iunit = 1, nunit
119 107444 : junit = junit + 1
120 476214 : r(1:3, iunit, imol) = particles%els(junit)%r(1:3)
121 : END DO
122 : END DO
123 :
124 : ! now let's find the LJ energy between all the oxygens and
125 : ! the charge interactions
126 8650 : lall = .TRUE.
127 8650 : jstart = 1
128 8650 : IF (PRESENT(molecule_number)) THEN
129 1808 : lall = .FALSE.
130 1808 : nstart = molecule_number
131 1808 : nend = molecule_number
132 : ELSE
133 14584 : nstart = 1
134 14584 : nend = SUM(nchains(:))
135 : END IF
136 30992 : DO imol = nstart, nend
137 23626 : IF (lall) jstart = imol + 1
138 23626 : nunits_i = nunits(mol_type(imol))
139 159526 : DO jmol = jstart, SUM(nchains(:))
140 93258 : IF (imol == jmol) CYCLE
141 91452 : nunits_j = nunits(mol_type(jmol))
142 :
143 350352 : DO iunit = 1, nunits_i
144 754516 : DO junit = 1, nunits_j
145 : ! find the minimum image distance
146 : RIJ(1) = r(1, iunit, imol) - r(1, junit, jmol) - &
147 : box_length(1)*ANINT( &
148 425984 : (r(1, iunit, imol) - r(1, junit, jmol))/box_length(1))
149 : RIJ(2) = r(2, iunit, imol) - r(2, junit, jmol) - &
150 : box_length(2)*ANINT( &
151 425984 : (r(2, iunit, imol) - r(2, junit, jmol))/box_length(2))
152 : RIJ(3) = r(3, iunit, imol) - r(3, junit, jmol) - &
153 : box_length(3)*ANINT( &
154 425984 : (r(3, iunit, imol) - r(3, junit, jmol))/box_length(3))
155 :
156 425984 : dist = RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2
157 :
158 662542 : IF (dist < rmin) THEN
159 1284 : loverlap = .TRUE.
160 1284 : DEALLOCATE (r)
161 :
162 1284 : CALL timestop(handle)
163 : RETURN
164 : END IF
165 :
166 : END DO
167 : END DO
168 : END DO
169 : END DO
170 :
171 7366 : DEALLOCATE (r)
172 :
173 : ! end the timing
174 7366 : CALL timestop(handle)
175 :
176 8650 : END SUBROUTINE check_for_overlap
177 :
178 : ! **************************************************************************************************
179 : !> \brief calculates the center of mass of a given molecule
180 : !> \param coordinates the coordinates of the atoms in the molecule
181 : !> \param natom the number of atoms in the molecule
182 : !> \param center_of_mass the coordinates of the center of mass
183 : !> \param mass the mass of the atoms in the molecule
184 : !>
185 : !> Designed for parallel use.
186 : !> \author MJM
187 : ! **************************************************************************************************
188 3432 : SUBROUTINE get_center_of_mass(coordinates, natom, center_of_mass, &
189 1144 : mass)
190 :
191 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: coordinates
192 : INTEGER, INTENT(IN) :: natom
193 : REAL(KIND=dp), DIMENSION(1:3), INTENT(OUT) :: center_of_mass
194 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: mass
195 :
196 : CHARACTER(len=*), PARAMETER :: routineN = 'get_center_of_mass'
197 :
198 : INTEGER :: handle, i, iatom
199 : REAL(KIND=dp) :: total_mass
200 :
201 : ! begin the timing of the subroutine
202 :
203 1144 : CALL timeset(routineN, handle)
204 :
205 3544 : total_mass = SUM(mass(1:natom))
206 1144 : center_of_mass(:) = 0.0E0_dp
207 :
208 3544 : DO iatom = 1, natom
209 10744 : DO i = 1, 3
210 : center_of_mass(i) = center_of_mass(i) + &
211 9600 : mass(iatom)*coordinates(i, iatom)
212 : END DO
213 : END DO
214 :
215 4576 : center_of_mass(1:3) = center_of_mass(1:3)/total_mass
216 :
217 : ! end the timing
218 1144 : CALL timestop(handle)
219 :
220 1144 : END SUBROUTINE get_center_of_mass
221 :
222 : ! **************************************************************************************************
223 : !> \brief folds all the coordinates into the center simulation box using
224 : !> a center of mass cutoff
225 : !> \param coordinates the coordinates of the atoms in the system
226 : !> \param nchains_tot the total number of molecules in the box
227 : !> \param mol_type an array that indicates the type of every molecule in the box
228 : !> \param mass the mass of every atom for all molecule types
229 : !> \param nunits the number of interaction sites for each molecule type
230 : !> \param box_length an array for the lengths of the simulation box sides
231 : !>
232 : !> Designed for parallel use.
233 : !> \author MJM
234 : ! **************************************************************************************************
235 0 : SUBROUTINE mc_coordinate_fold(coordinates, nchains_tot, mol_type, mass, nunits, box_length)
236 :
237 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: coordinates
238 : INTEGER, INTENT(IN) :: nchains_tot
239 : INTEGER, DIMENSION(:), INTENT(IN) :: mol_type
240 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: mass
241 : INTEGER, DIMENSION(:), INTENT(IN) :: nunits
242 : REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: box_length
243 :
244 : CHARACTER(len=*), PARAMETER :: routineN = 'mc_coordinate_fold'
245 :
246 : INTEGER :: end_atom, handle, iatom, imolecule, &
247 : jatom, molecule_type, natoms, &
248 : start_atom
249 : REAL(KIND=dp), DIMENSION(1:3) :: center_of_mass
250 :
251 : ! begin the timing of the subroutine
252 :
253 0 : CALL timeset(routineN, handle)
254 :
255 : ! loop over all molecules
256 0 : end_atom = 0
257 0 : DO imolecule = 1, nchains_tot
258 0 : molecule_type = mol_type(imolecule)
259 0 : natoms = nunits(molecule_type)
260 0 : start_atom = end_atom + 1
261 0 : end_atom = start_atom + natoms - 1
262 : CALL get_center_of_mass(coordinates(:, start_atom:end_atom), &
263 0 : natoms, center_of_mass(:), mass(:, molecule_type))
264 0 : DO iatom = 1, natoms
265 0 : jatom = iatom + start_atom - 1
266 : coordinates(1, jatom) = coordinates(1, jatom) - &
267 0 : box_length(1)*FLOOR(center_of_mass(1)/box_length(1))
268 : coordinates(2, jatom) = coordinates(2, jatom) - &
269 0 : box_length(2)*FLOOR(center_of_mass(2)/box_length(2))
270 : coordinates(3, jatom) = coordinates(3, jatom) - &
271 0 : box_length(3)*FLOOR(center_of_mass(3)/box_length(2))
272 : END DO
273 :
274 : END DO
275 :
276 : ! end the timing
277 0 : CALL timestop(handle)
278 :
279 0 : END SUBROUTINE mc_coordinate_fold
280 :
281 : ! **************************************************************************************************
282 : !> \brief takes the last molecule in a force environment and moves it around
283 : !> to different center of mass positions and orientations, selecting one
284 : !> based on the rosenbluth weight
285 : !> \param force_env the force environment containing the coordinates
286 : !> \param BETA the value of 1/kT for this simulations, in a.u.
287 : !> \param max_val ...
288 : !> \param min_val ...
289 : !> \param exp_max_val ...
290 : !> \param exp_min_val ...
291 : !> \param nswapmoves the number of desired trial configurations
292 : !> \param rosenbluth_weight the Rosenbluth weight for this set of configs
293 : !> \param start_atom the atom number that the molecule to be swapped starts on
294 : !> \param natoms_tot the total number of interaction sites in the box
295 : !> \param nunits the number of interaction sites for every molecule_type
296 : !> \param nunits_mol ...
297 : !> \param mass the mass for every interaction site of every molecule type
298 : !> \param loverlap the flag that indicates if all of the configs have an
299 : !> atomic overlap
300 : !> \param choosen_energy the energy of the chosen config
301 : !> \param old_energy the energy that we subtract from all of the trial
302 : !> energies to prevent numerical overflows
303 : !> \param ionode indicates if we're on the main CPU
304 : !> \param lremove is this the Rosenbluth weight for a removal box?
305 : !> \param mol_type an array that contains the molecule type for every atom in the box
306 : !> \param nchains the number of molecules of each type in this box
307 : !> \param source the MPI source value, for broadcasts
308 : !> \param group the MPI group value, for broadcasts
309 : !> \param rng_stream the random number stream that we draw from
310 : !> \param avbmc_atom ...
311 : !> \param rmin ...
312 : !> \param rmax ...
313 : !> \param move_type ...
314 : !> \param target_atom ...
315 : !> \par Optional Avbmc Flags
316 : !> - avbmc_atom: the atom number that serves for the target atom in each
317 : !> molecule (1 is the first atom in the molecule, etc.)
318 : !> - rmin: the minimum AVBMC radius for the shell around the target
319 : !> - rmax: the maximum AVBMC radius for the shell around the target
320 : !> - move_type: generate configs in the "in" or "out" volume
321 : !> - target_atom: the number of the avbmc atom in the target molecule
322 : !> \par
323 : !> Suitable for parallel.
324 : !> \author MJM
325 : ! **************************************************************************************************
326 44 : SUBROUTINE generate_cbmc_swap_config(force_env, BETA, max_val, min_val, exp_max_val, &
327 88 : exp_min_val, nswapmoves, rosenbluth_weight, start_atom, natoms_tot, nunits, nunits_mol, &
328 44 : mass, loverlap, choosen_energy, old_energy, ionode, lremove, mol_type, nchains, source, &
329 : group, rng_stream, avbmc_atom, rmin, rmax, move_type, target_atom)
330 :
331 : TYPE(force_env_type), POINTER :: force_env
332 : REAL(KIND=dp), INTENT(IN) :: BETA, max_val, min_val, exp_max_val, &
333 : exp_min_val
334 : INTEGER, INTENT(IN) :: nswapmoves
335 : REAL(KIND=dp), INTENT(OUT) :: rosenbluth_weight
336 : INTEGER, INTENT(IN) :: start_atom, natoms_tot
337 : INTEGER, DIMENSION(:), INTENT(IN) :: nunits
338 : INTEGER, INTENT(IN) :: nunits_mol
339 : REAL(dp), DIMENSION(1:nunits_mol), INTENT(IN) :: mass
340 : LOGICAL, INTENT(OUT) :: loverlap
341 : REAL(KIND=dp), INTENT(OUT) :: choosen_energy
342 : REAL(KIND=dp), INTENT(IN) :: old_energy
343 : LOGICAL, INTENT(IN) :: ionode, lremove
344 : INTEGER, DIMENSION(:), INTENT(IN) :: mol_type, nchains
345 : INTEGER, INTENT(IN) :: source
346 : TYPE(mp_comm_type) :: group
347 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
348 : INTEGER, INTENT(IN), OPTIONAL :: avbmc_atom
349 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rmin, rmax
350 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: move_type
351 : INTEGER, INTENT(IN), OPTIONAL :: target_atom
352 :
353 : CHARACTER(len=*), PARAMETER :: routineN = 'generate_cbmc_swap_config'
354 :
355 : INTEGER :: atom_number, choosen, end_atom, handle, &
356 : i, iatom, imolecule, imove, &
357 : molecule_number
358 : LOGICAL :: all_overlaps
359 44 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: loverlap_array
360 : REAL(KIND=dp) :: bias_energy, exponent, rand, &
361 : total_running_weight
362 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: boltz_weights, trial_energy
363 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: r_old
364 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r
365 : REAL(KIND=dp), DIMENSION(1:3) :: abc, center_of_mass, diff, r_insert
366 : TYPE(cell_type), POINTER :: cell
367 : TYPE(cp_subsys_type), POINTER :: oldsys
368 : TYPE(particle_list_type), POINTER :: particles
369 :
370 : ! begin the timing of the subroutine
371 :
372 44 : CALL timeset(routineN, handle)
373 :
374 44 : NULLIFY (oldsys)
375 : ! get the particle coordinates and the cell length
376 44 : CALL force_env_get(force_env, cell=cell, subsys=oldsys)
377 44 : CALL get_cell(cell, abc=abc)
378 44 : CALL cp_subsys_get(oldsys, particles=particles)
379 :
380 : ! do some checking to make sure we have all the data we need
381 44 : IF (PRESENT(avbmc_atom)) THEN
382 : IF (.NOT. PRESENT(rmin) .OR. .NOT. PRESENT(rmax) .OR. &
383 0 : .NOT. PRESENT(move_type) .OR. .NOT. PRESENT(target_atom)) THEN
384 0 : CPABORT('AVBMC swap move is missing information!')
385 : END IF
386 : END IF
387 :
388 132 : ALLOCATE (r_old(1:3, 1:natoms_tot))
389 176 : ALLOCATE (r(1:3, 1:natoms_tot, 1:nswapmoves))
390 132 : ALLOCATE (trial_energy(1:nswapmoves))
391 132 : ALLOCATE (boltz_weights(1:nswapmoves))
392 132 : ALLOCATE (loverlap_array(1:nswapmoves))
393 :
394 : ! initialize the arrays that need it
395 220 : loverlap_array(:) = .FALSE.
396 44 : loverlap = .FALSE.
397 220 : boltz_weights(:) = 0.0_dp
398 220 : trial_energy(:) = 0.0_dp
399 18492 : r(:, :, :) = 0.0_dp
400 44 : choosen_energy = 0.0_dp
401 44 : rosenbluth_weight = 0.0_dp
402 :
403 : ! save the positions of the molecules
404 220 : DO imove = 1, nswapmoves
405 4788 : DO iatom = 1, natoms_tot
406 18448 : r(1:3, iatom, imove) = particles%els(iatom)%r(1:3)
407 : END DO
408 : END DO
409 :
410 : ! save the remove coordinates
411 1186 : DO iatom = 1, natoms_tot
412 4612 : r_old(1:3, iatom) = r(1:3, iatom, 1)
413 : END DO
414 :
415 : ! figure out the numbers of the first and last atoms in the molecule
416 44 : end_atom = start_atom + nunits_mol - 1
417 : ! figure out which molecule number we're on
418 44 : molecule_number = 0
419 44 : atom_number = 1
420 456 : DO imolecule = 1, SUM(nchains(:))
421 368 : IF (atom_number == start_atom) THEN
422 44 : molecule_number = imolecule
423 44 : EXIT
424 : END IF
425 324 : atom_number = atom_number + nunits(mol_type(imolecule))
426 : END DO
427 44 : IF (molecule_number == 0) CALL cp_abort(__LOCATION__, &
428 0 : 'CBMC swap move cannot find which molecule number it needs')
429 :
430 44 : IF (lremove) THEN
431 : CALL check_for_overlap(force_env, nchains, nunits, loverlap_array(1), &
432 22 : mol_type)
433 22 : CALL group%bcast(loverlap_array(1), source)
434 :
435 22 : IF (loverlap_array(1)) THEN
436 0 : IF (ionode) THEN
437 0 : WRITE (*, *) start_atom, end_atom, natoms_tot
438 0 : DO iatom = 1, natoms_tot
439 0 : WRITE (*, *) r(1:3, iatom, 1)
440 : END DO
441 : END IF
442 0 : CPABORT('CBMC swap move found an overlap in the old config')
443 : END IF
444 : END IF
445 :
446 220 : DO imove = 1, nswapmoves
447 :
448 : ! drop into serial
449 176 : IF (ionode) THEN
450 :
451 88 : IF (PRESENT(avbmc_atom)) THEN
452 : ! find an AVBMC insertion point
453 : CALL generate_avbmc_insertion(rmin, rmax, &
454 : r_old(1:3, target_atom), &
455 0 : move_type, r_insert(:), abc(:), rng_stream)
456 :
457 0 : DO i = 1, 3
458 0 : diff(i) = r_insert(i) - r_old(i, start_atom + avbmc_atom - 1)
459 : END DO
460 :
461 : ELSE
462 : ! find a new insertion point somewhere in the box
463 352 : DO i = 1, 3
464 264 : rand = rng_stream%next()
465 352 : r_insert(i) = rand*abc(i)
466 : END DO
467 :
468 : ! find the center of mass of the insertion molecule
469 : CALL get_center_of_mass(r(:, start_atom:end_atom, imove), nunits_mol, &
470 88 : center_of_mass(:), mass(:))
471 :
472 : ! move the molecule to the insertion point
473 :
474 352 : DO i = 1, 3
475 352 : diff(i) = r_insert(i) - center_of_mass(i)
476 : END DO
477 :
478 : END IF
479 :
480 256 : DO iatom = start_atom, end_atom
481 760 : r(1:3, iatom, imove) = r(1:3, iatom, imove) + diff(1:3)
482 : END DO
483 :
484 : ! rotate the molecule...this routine is only made for serial use
485 : CALL rotate_molecule(r(:, start_atom:end_atom, imove), mass(:), &
486 88 : nunits_mol, rng_stream)
487 :
488 88 : IF (imove == 1 .AND. lremove) THEN
489 293 : DO iatom = 1, natoms_tot
490 1139 : r(1:3, iatom, 1) = r_old(1:3, iatom)
491 : END DO
492 : END IF
493 :
494 : END IF
495 :
496 176 : CALL group%bcast(r(:, :, imove), source)
497 :
498 : ! calculate the energy and boltzman weight of the config
499 512 : DO iatom = start_atom, end_atom
500 1520 : particles%els(iatom)%r(1:3) = r(1:3, iatom, imove)
501 : END DO
502 :
503 : CALL check_for_overlap(force_env, nchains, nunits, loverlap_array(imove), &
504 176 : mol_type, molecule_number=molecule_number)
505 176 : IF (loverlap_array(imove)) THEN
506 6 : boltz_weights(imove) = 0.0_dp
507 6 : CYCLE
508 : END IF
509 :
510 170 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
511 : CALL force_env_get(force_env, &
512 170 : potential_energy=bias_energy)
513 :
514 170 : trial_energy(imove) = (bias_energy - old_energy)
515 170 : exponent = -BETA*trial_energy(imove)
516 :
517 214 : IF (exponent .GT. exp_max_val) THEN
518 0 : boltz_weights(imove) = max_val
519 170 : ELSEIF (exponent .LT. exp_min_val) THEN
520 2 : boltz_weights(imove) = min_val
521 : ELSE
522 168 : boltz_weights(imove) = EXP(exponent)
523 : END IF
524 :
525 : END DO
526 :
527 : ! now we need to pick a configuration based on the Rosenbluth weight,
528 : ! which is just the sum of the Boltzmann weights
529 220 : rosenbluth_weight = SUM(boltz_weights(:))
530 44 : IF (rosenbluth_weight .EQ. 0.0_dp .AND. lremove) THEN
531 : ! should never have 0.0 for an old weight...causes a divide by zero
532 : ! in the acceptance rule
533 0 : IF (ionode) THEN
534 0 : WRITE (*, *) boltz_weights(1:nswapmoves)
535 0 : WRITE (*, *) start_atom, end_atom, lremove
536 0 : WRITE (*, *) loverlap_array(1:nswapmoves)
537 0 : WRITE (*, *) natoms_tot
538 0 : WRITE (*, *)
539 0 : DO iatom = 1, natoms_tot
540 0 : WRITE (*, *) r(1:3, iatom, 1)*angstrom
541 : END DO
542 : END IF
543 0 : CPABORT('CBMC swap move found a bad old weight')
544 : END IF
545 44 : all_overlaps = .TRUE.
546 44 : total_running_weight = 0.0E0_dp
547 44 : choosen = 0
548 44 : IF (ionode) THEN
549 22 : rand = rng_stream%next()
550 : ! CALL random_number(rand)
551 : END IF
552 44 : CALL group%bcast(rand, source)
553 98 : DO imove = 1, nswapmoves
554 98 : IF (loverlap_array(imove)) CYCLE
555 94 : all_overlaps = .FALSE.
556 94 : total_running_weight = total_running_weight + boltz_weights(imove)
557 94 : IF (total_running_weight .GE. rand*rosenbluth_weight) THEN
558 : choosen = imove
559 : EXIT
560 : END IF
561 : END DO
562 :
563 0 : IF (all_overlaps) THEN
564 0 : loverlap = .TRUE.
565 :
566 : ! if this is an old configuration, we always choose the first one...
567 : ! this should never be the case, but I'm testing something
568 0 : IF (lremove) THEN
569 0 : IF (ionode) THEN
570 0 : WRITE (*, *) boltz_weights(1:nswapmoves)
571 0 : WRITE (*, *) start_atom, end_atom, lremove
572 0 : WRITE (*, *) loverlap_array(1:nswapmoves)
573 0 : DO iatom = 1, natoms_tot
574 0 : WRITE (*, *) r(1:3, iatom, 1)
575 : END DO
576 : END IF
577 0 : CPABORT('CBMC swap move found all overlaps for the remove config')
578 : END IF
579 :
580 0 : DEALLOCATE (r_old)
581 0 : DEALLOCATE (r)
582 0 : DEALLOCATE (trial_energy)
583 0 : DEALLOCATE (boltz_weights)
584 0 : DEALLOCATE (loverlap_array)
585 0 : CALL timestop(handle)
586 0 : RETURN
587 : END IF
588 :
589 : ! make sure a configuration was chosen
590 44 : IF (choosen == 0) &
591 0 : CPABORT('CBMC swap move failed to select config')
592 :
593 : ! if this is an old configuration, we always choose the first one
594 44 : IF (lremove) choosen = 1
595 :
596 : ! set the energy for the configuration
597 44 : choosen_energy = trial_energy(choosen)
598 :
599 : ! copy the coordinates to the force environment
600 1186 : DO iatom = 1, natoms_tot
601 4612 : particles%els(iatom)%r(1:3) = r(1:3, iatom, choosen)
602 : END DO
603 :
604 44 : DEALLOCATE (r_old)
605 44 : DEALLOCATE (r)
606 44 : DEALLOCATE (trial_energy)
607 44 : DEALLOCATE (boltz_weights)
608 44 : DEALLOCATE (loverlap_array)
609 :
610 : ! end the timing
611 44 : CALL timestop(handle)
612 :
613 44 : END SUBROUTINE generate_cbmc_swap_config
614 :
615 : ! **************************************************************************************************
616 : !> \brief rotates a molecule randomly around the center of mass,
617 : !> sequentially in x, y, and z directions
618 : !> \param r the coordinates of the molecule to rotate
619 : !> \param mass the mass of all the atoms in the molecule
620 : !> \param natoms the number of atoms in the molecule
621 : !> \param rng_stream the stream we pull random numbers from
622 : !>
623 : !> Use only in serial.
624 : !> \author MJM
625 : ! **************************************************************************************************
626 98 : SUBROUTINE rotate_molecule(r, mass, natoms, rng_stream)
627 :
628 : INTEGER, INTENT(IN) :: natoms
629 : REAL(KIND=dp), DIMENSION(1:natoms), INTENT(IN) :: mass
630 : REAL(KIND=dp), DIMENSION(1:3, 1:natoms), &
631 : INTENT(INOUT) :: r
632 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
633 :
634 : CHARACTER(len=*), PARAMETER :: routineN = 'rotate_molecule'
635 :
636 : INTEGER :: handle, iunit
637 : REAL(KIND=dp) :: cosdg, dgamma, rand, rx, rxnew, ry, &
638 : rynew, rz, rznew, sindg
639 : REAL(KIND=dp), DIMENSION(1:3) :: center_of_mass
640 :
641 : ! begin the timing of the subroutine
642 :
643 98 : CALL timeset(routineN, handle)
644 :
645 : ! find the center of mass of the molecule
646 98 : CALL get_center_of_mass(r(:, :), natoms, center_of_mass(:), mass(:))
647 :
648 : ! call a random number to figure out how far we're moving
649 98 : rand = rng_stream%next()
650 98 : dgamma = pi*(rand - 0.5E0_dp)*2.0E0_dp
651 :
652 : ! *** set up the rotation matrix ***
653 :
654 98 : cosdg = COS(dgamma)
655 98 : sindg = SIN(dgamma)
656 :
657 : ! *** ROTATE UNITS OF I AROUND X-AXIS ***
658 :
659 296 : DO iunit = 1, natoms
660 198 : ry = r(2, iunit) - center_of_mass(2)
661 198 : rz = r(3, iunit) - center_of_mass(3)
662 198 : rynew = cosdg*ry + sindg*rz
663 198 : rznew = cosdg*rz - sindg*ry
664 :
665 198 : r(2, iunit) = rynew + center_of_mass(2)
666 296 : r(3, iunit) = rznew + center_of_mass(3)
667 :
668 : END DO
669 :
670 : ! *** ROTATE UNITS OF I AROUND y-AXIS ***
671 :
672 296 : DO iunit = 1, natoms
673 198 : rx = r(1, iunit) - center_of_mass(1)
674 198 : rz = r(3, iunit) - center_of_mass(3)
675 198 : rxnew = cosdg*rx + sindg*rz
676 198 : rznew = cosdg*rz - sindg*rx
677 :
678 198 : r(1, iunit) = rxnew + center_of_mass(1)
679 296 : r(3, iunit) = rznew + center_of_mass(3)
680 :
681 : END DO
682 :
683 : ! *** ROTATE UNITS OF I AROUND z-AXIS ***
684 :
685 296 : DO iunit = 1, natoms
686 198 : rx = r(1, iunit) - center_of_mass(1)
687 198 : ry = r(2, iunit) - center_of_mass(2)
688 198 : rxnew = cosdg*rx + sindg*ry
689 198 : rynew = cosdg*ry - sindg*rx
690 :
691 198 : r(1, iunit) = rxnew + center_of_mass(1)
692 296 : r(2, iunit) = rynew + center_of_mass(2)
693 :
694 : END DO
695 :
696 : ! end the timing
697 98 : CALL timestop(handle)
698 :
699 98 : END SUBROUTINE rotate_molecule
700 :
701 : ! **************************************************************************************************
702 : !> \brief selects a molecule at random to perform a MC move on...you can specify
703 : !> the box the molecule should be in, its type, both, or neither
704 : !> \param mc_molecule_info the structure that contains some global molecule data
705 : !> \param start_atom the number of the first atom in the chosen molecule in relation
706 : !> to the force_env it's in
707 : !> \param box_number the box the chosen molecule is in
708 : !> \param molecule_type the type of molecule the chosen molecule is
709 : !> \param rng_stream the stream we pull random numbers from
710 : !> \param box if present, tells the routine which box to grab a molecule from
711 : !> \param molecule_type_old if present, tells the routine which molecule type to select from
712 : !> \author MJM
713 : ! **************************************************************************************************
714 816 : SUBROUTINE find_mc_test_molecule(mc_molecule_info, start_atom, &
715 : box_number, molecule_type, rng_stream, box, molecule_type_old)
716 :
717 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
718 : INTEGER, INTENT(OUT) :: start_atom, box_number, molecule_type
719 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
720 : INTEGER, INTENT(IN), OPTIONAL :: box, molecule_type_old
721 :
722 : CHARACTER(LEN=*), PARAMETER :: routineN = 'find_mc_test_molecule'
723 :
724 : INTEGER :: handle, ibox, imol_type, imolecule, &
725 : jbox, molecule_number, nchains_tot, &
726 : start_mol
727 816 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits
728 816 : INTEGER, DIMENSION(:, :), POINTER :: nchains
729 : REAL(KIND=dp) :: rand
730 :
731 : ! begin the timing of the subroutine
732 :
733 816 : CALL timeset(routineN, handle)
734 :
735 816 : NULLIFY (nunits, mol_type, nchains)
736 : CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, nunits=nunits, &
737 816 : mol_type=mol_type)
738 :
739 : ! initialize the outgoing variables
740 816 : start_atom = 0
741 816 : box_number = 0
742 816 : molecule_type = 0
743 :
744 816 : IF (PRESENT(box) .AND. PRESENT(molecule_type_old)) THEN
745 : ! only need to find the atom number the molecule starts on
746 253 : rand = rng_stream%next()
747 253 : molecule_number = CEILING(rand*REAL(nchains(molecule_type_old, box), KIND=dp))
748 :
749 253 : start_mol = 1
750 253 : DO jbox = 1, box - 1
751 253 : start_mol = start_mol + SUM(nchains(:, jbox))
752 : END DO
753 :
754 : ! adjust to take into account molecules of other types in the box
755 253 : DO imol_type = 1, molecule_type_old - 1
756 253 : molecule_number = molecule_number + nchains(imol_type, box)
757 : END DO
758 :
759 253 : start_atom = 1
760 981 : DO imolecule = 1, molecule_number - 1
761 981 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
762 : END DO
763 :
764 563 : ELSEIF (PRESENT(box)) THEN
765 : ! any molecule in box...need to find molecule type and start atom
766 0 : rand = rng_stream%next()
767 0 : molecule_number = CEILING(rand*REAL(SUM(nchains(:, box)), KIND=dp))
768 :
769 0 : start_mol = 1
770 0 : DO jbox = 1, box - 1
771 0 : start_mol = start_mol + SUM(nchains(:, jbox))
772 : END DO
773 :
774 0 : molecule_type = mol_type(start_mol + molecule_number - 1)
775 :
776 : ! now the starting atom
777 0 : start_atom = 1
778 0 : DO imolecule = 1, molecule_number - 1
779 0 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
780 : END DO
781 :
782 563 : ELSEIF (PRESENT(molecule_type_old)) THEN
783 : ! any molecule of type molecule_type_old...need to find box number and start atom
784 563 : rand = rng_stream%next()
785 1162 : molecule_number = CEILING(rand*REAL(SUM(nchains(molecule_type_old, :)), KIND=dp))
786 :
787 : ! find which box it's in
788 563 : nchains_tot = 0
789 578 : DO ibox = 1, SIZE(nchains(molecule_type_old, :))
790 578 : IF (molecule_number .LE. nchains(molecule_type_old, ibox)) THEN
791 563 : box_number = ibox
792 563 : EXIT
793 : END IF
794 15 : molecule_number = molecule_number - nchains(molecule_type_old, ibox)
795 : END DO
796 :
797 563 : start_mol = 1
798 578 : DO jbox = 1, box_number - 1
799 608 : start_mol = start_mol + SUM(nchains(:, jbox))
800 : END DO
801 :
802 : ! now find the starting atom number
803 713 : DO imol_type = 1, molecule_type_old - 1
804 713 : molecule_number = molecule_number + nchains(imol_type, box_number)
805 : END DO
806 563 : start_atom = 1
807 3237 : DO imolecule = 1, molecule_number - 1
808 3237 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
809 : END DO
810 :
811 : ELSE
812 : ! no restrictions...need to find all pieces of data
813 0 : nchains_tot = 0
814 0 : DO ibox = 1, SIZE(nchains(1, :))
815 0 : nchains_tot = nchains_tot + SUM(nchains(:, ibox))
816 : END DO
817 0 : rand = rng_stream%next()
818 0 : molecule_number = CEILING(rand*REAL(nchains_tot, KIND=dp))
819 :
820 0 : molecule_type = mol_type(molecule_number)
821 :
822 : ! now which box it's in
823 0 : DO ibox = 1, SUM(nchains(1, :))
824 0 : IF (molecule_number .LE. SUM(nchains(:, ibox))) THEN
825 0 : box_number = ibox
826 0 : EXIT
827 : END IF
828 0 : molecule_number = molecule_number - SUM(nchains(:, ibox))
829 : END DO
830 :
831 : ! now find the starting atom number
832 0 : start_mol = 1
833 0 : DO jbox = 1, box_number - 1
834 0 : start_mol = start_mol + SUM(nchains(:, jbox))
835 : END DO
836 0 : start_atom = 1
837 0 : DO imolecule = 1, molecule_number - 1
838 0 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
839 : END DO
840 :
841 : END IF
842 :
843 : ! make sure things are good
844 816 : IF (PRESENT(box)) box_number = box
845 816 : IF (PRESENT(molecule_type_old)) molecule_type = molecule_type_old
846 :
847 816 : CPASSERT(start_atom /= 0)
848 816 : CPASSERT(box_number /= 0)
849 816 : CPASSERT(molecule_type /= 0)
850 :
851 : ! end the timing
852 816 : CALL timestop(handle)
853 :
854 816 : END SUBROUTINE find_mc_test_molecule
855 :
856 : ! **************************************************************************************************
857 : !> \brief generates an array that tells us which sides of the simulation
858 : !> cell we can increase or decrease using a discrete volume move
859 : !> \param cell the lengths of the sides of the cell
860 : !> \param discrete_array the array that indicates which sides we can move
861 : !> \param step_size the size of the discrete volume move
862 : !>
863 : !> Suitable for parallel.
864 : !> \author MJM
865 : ! **************************************************************************************************
866 0 : SUBROUTINE create_discrete_array(cell, discrete_array, step_size)
867 :
868 : ! 1 is for increase, 2 is for decrease
869 : ! 1 is for "yes, we can do the move", 0 is for no
870 :
871 : REAL(dp), DIMENSION(1:3), INTENT(IN) :: cell
872 : INTEGER, DIMENSION(1:3, 1:2), INTENT(OUT) :: discrete_array
873 : REAL(dp), INTENT(IN) :: step_size
874 :
875 : INTEGER :: iside
876 : REAL(dp) :: high_value, length1, length2, low_value
877 :
878 0 : discrete_array(:, :) = 0
879 :
880 0 : length1 = ABS(cell(1) - cell(2))
881 0 : length2 = ABS(cell(2) - cell(3))
882 :
883 : ! now let's figure out all the different cases
884 0 : IF (length1 .LT. 0.01_dp*step_size .AND. &
885 : length2 .LT. 0.01_dp*step_size) THEN
886 : ! all sides are equal, so we can move up or down
887 0 : discrete_array(1:3, 1) = 1
888 0 : discrete_array(1:3, 2) = 1
889 : ELSE
890 :
891 : ! find the low value and the high value
892 0 : high_value = -1.0_dp
893 0 : low_value = cell(1)*cell(2)*cell(3)
894 0 : DO iside = 1, 3
895 0 : IF (cell(iside) .LT. low_value) low_value = cell(iside)
896 0 : IF (cell(iside) .GT. high_value) high_value = cell(iside)
897 : END DO
898 0 : DO iside = 1, 3
899 : ! now we see if the value is a high value or a low value...it can only be
900 : ! one of the two
901 0 : IF (ABS(cell(iside) - low_value) .LT. 0.01_dp*step_size) THEN
902 : ! low value, we can only increase the cell size
903 0 : discrete_array(iside, 1) = 1
904 0 : discrete_array(iside, 2) = 0
905 : ELSE
906 : ! high value, we can only decrease the cell size
907 0 : discrete_array(iside, 1) = 0
908 0 : discrete_array(iside, 2) = 1
909 : END IF
910 : END DO
911 : END IF
912 :
913 0 : END SUBROUTINE create_discrete_array
914 :
915 : ! **************************************************************************************************
916 : !> \brief generates an insertion point in either the "in" or the "out" volume
917 : !> of a target atom, where the "in" volume is a shell with inner radius
918 : !> rmin and outer radius rmax
919 : !> \param rmin the minimum AVBMC radius for the shell around the target
920 : !> \param rmax the maximum AVBMC radius for the shell around the target
921 : !> \param r_target the coordinates of the target atom
922 : !> \param move_type generate configs in the "in" or "out" volume
923 : !> \param r_insert the output insertion site
924 : !> \param abc the lengths of the sides of the simulation box
925 : !> \param rng_stream the random number stream that we draw from
926 : !>
927 : !> Use only in serial.
928 : !> \author MJM
929 : ! **************************************************************************************************
930 0 : SUBROUTINE generate_avbmc_insertion(rmin, rmax, r_target, &
931 : move_type, r_insert, abc, rng_stream)
932 :
933 : REAL(KIND=dp), INTENT(IN) :: rmin, rmax
934 : REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: r_target
935 : CHARACTER(LEN=*), INTENT(IN) :: move_type
936 : REAL(KIND=dp), DIMENSION(1:3), INTENT(OUT) :: r_insert
937 : REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: abc
938 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
939 :
940 : INTEGER :: i
941 : REAL(dp) :: dist, eta_1, eta_2, eta_sq, rand
942 : REAL(dp), DIMENSION(1:3) :: RIJ
943 :
944 0 : r_insert(1:3) = 0.0_dp
945 :
946 0 : IF (move_type == 'in') THEN
947 : ! generate a random unit vector, from Allen and Tildesley
948 : DO
949 0 : eta_1 = rng_stream%next()
950 0 : eta_2 = rng_stream%next()
951 0 : eta_sq = eta_1**2 + eta_2**2
952 0 : IF (eta_sq .LT. 1.0_dp) THEN
953 0 : r_insert(1) = 2.0_dp*eta_1*SQRT(1.0_dp - eta_sq)
954 0 : r_insert(2) = 2.0_dp*eta_2*SQRT(1.0_dp - eta_sq)
955 0 : r_insert(3) = 1.0_dp - 2.0_dp*eta_sq
956 : EXIT
957 : END IF
958 : END DO
959 :
960 : ! now scale that vector to be within the "in" region
961 0 : rand = rng_stream%next()
962 : r_insert(1:3) = r_insert(1:3)*(rand*(rmax**3 - rmin**3) + rmin**3)** &
963 0 : (1.0_dp/3.0_dp)
964 :
965 0 : r_insert(1:3) = r_target(1:3) + r_insert(1:3)
966 : ELSE
967 :
968 : ! find a new insertion point somewhere in the box
969 : DO
970 0 : DO i = 1, 3
971 0 : rand = rng_stream%next()
972 0 : r_insert(i) = rand*abc(i)
973 : END DO
974 :
975 : ! make sure it's not in the "in" region
976 : RIJ(1) = r_insert(1) - r_target(1) - abc(1)* &
977 0 : ANINT((r_insert(1) - r_target(1))/abc(1))
978 : RIJ(2) = r_insert(2) - r_target(2) - abc(2)* &
979 0 : ANINT((r_insert(2) - r_target(2))/abc(2))
980 : RIJ(3) = r_insert(3) - r_target(3) - abc(3)* &
981 0 : ANINT((r_insert(3) - r_target(3))/abc(3))
982 :
983 0 : dist = RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2
984 :
985 0 : IF (dist .LT. rmin**2 .OR. dist .GT. rmax**2) THEN
986 : EXIT
987 : END IF
988 :
989 : END DO
990 : END IF
991 :
992 0 : END SUBROUTINE generate_avbmc_insertion
993 :
994 : ! *****************************************************************************
995 : ! **************************************************************************************************
996 : !> \brief determine the number of cluster present in the given configuration
997 : !> based on the rclus value
998 : !> \param mc_par the mc parameters for the force env
999 : !> \param force_env the force environment containing the coordinates
1000 : !> \param cluster ...
1001 : !> \param nchains the number of molecules of each type in the box
1002 : !> \param nunits the number of interaction sites for each molecule
1003 : !> \param mol_type an array that indicates the type of each molecule
1004 : !> \param total_clus ...
1005 : !> \par
1006 : !> Original Multiparticle/Cluster Translation paper:
1007 : !> Orkoulas, Gerassimos, and Athanassios Z. Panagiotopoulos. Free energy and
1008 : !> phase equilibria for the restricted primitive model of ionic fluids from Monte
1009 : !> Carlo simulations. J. Chem. Phys. 1994,101.2,1452-1459.
1010 : !> \author Himanshu Goel
1011 : ! **************************************************************************************************
1012 :
1013 20 : SUBROUTINE cluster_search(mc_par, force_env, cluster, nchains, nunits, mol_type, total_clus)
1014 :
1015 : TYPE(mc_simpar_type), POINTER :: mc_par
1016 : TYPE(force_env_type), POINTER :: force_env
1017 : INTEGER, DIMENSION(:, :), INTENT(INOUT) :: cluster
1018 : INTEGER, DIMENSION(:), INTENT(IN) :: nchains, nunits, mol_type
1019 : INTEGER, INTENT(INOUT) :: total_clus
1020 :
1021 : CHARACTER(len=*), PARAMETER :: routineN = 'cluster_search'
1022 :
1023 : INTEGER :: counter, handle, imol, iunit, jmol, &
1024 : junit, nend, nstart, nunit, nunits_i, &
1025 : nunits_j
1026 20 : INTEGER, ALLOCATABLE, DIMENSION(:) :: clusmat, decision
1027 : LOGICAL :: lclus
1028 : REAL(KIND=dp) :: dx, dy, dz, rclus, rclussquare, rsquare
1029 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xcoord, ycoord, zcoord
1030 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r
1031 : REAL(KIND=dp), DIMENSION(1:3) :: abc
1032 : TYPE(cell_type), POINTER :: cell
1033 : TYPE(cp_subsys_type), POINTER :: oldsys
1034 : TYPE(particle_list_type), POINTER :: particles
1035 :
1036 : ! begin the timing of the subroutine
1037 :
1038 20 : CALL timeset(routineN, handle)
1039 :
1040 20 : NULLIFY (oldsys, particles)
1041 :
1042 : ! Getting Particle Coordinates
1043 :
1044 20 : CALL force_env_get(force_env, cell=cell, subsys=oldsys)
1045 20 : CALL get_cell(cell, abc=abc)
1046 20 : CALL cp_subsys_get(oldsys, particles=particles)
1047 20 : CALL get_mc_par(mc_par, rclus=rclus)
1048 :
1049 120 : ALLOCATE (r(1:3, 1:MAXVAL(nunits), 1:SUM(nchains)))
1050 :
1051 : ! Arranging particles coordinates into an easier matrix
1052 40 : nend = SUM(nchains(:))
1053 : junit = 0
1054 120 : DO imol = 1, nend
1055 100 : nunit = nunits(mol_type(imol))
1056 320 : DO iunit = 1, nunit
1057 200 : junit = junit + 1
1058 900 : r(1:3, iunit, imol) = particles%els(junit)%r(1:3)
1059 : END DO
1060 : END DO
1061 :
1062 20 : counter = 0
1063 :
1064 : ! Allocating the size of matrix and decision matrix
1065 100 : ALLOCATE (clusmat(nend), decision(nend))
1066 :
1067 : !Initialize
1068 120 : DO imol = 1, nend
1069 100 : decision(imol) = 0
1070 120 : clusmat(imol) = 0
1071 : END DO
1072 :
1073 20 : rclussquare = rclus*rclus
1074 : ! Starting the cluster count loop
1075 720 : DO WHILE (SUM(decision) .LT. nend)
1076 300 : DO nstart = 1, nend
1077 300 : IF (clusmat(nstart) .EQ. 0) THEN
1078 100 : counter = counter + 1
1079 100 : clusmat(nstart) = counter
1080 100 : EXIT
1081 : END IF
1082 : END DO
1083 :
1084 100 : lclus = .TRUE.
1085 320 : DO WHILE (lclus .EQV. .TRUE.)
1086 900 : DO imol = 1, nend
1087 800 : nunits_i = nunits(mol_type(imol))
1088 : ! Allocating the xcoord,ycoord,zcoord based upon the size of molecule nunits
1089 800 : lclus = .FALSE.
1090 900 : IF (clusmat(imol) .EQ. counter .AND. decision(imol) .EQ. 0) THEN
1091 700 : ALLOCATE (xcoord(nunits_i), ycoord(nunits_i), zcoord(nunits_i))
1092 100 : decision(imol) = 1
1093 100 : lclus = .TRUE.
1094 300 : DO iunit = 1, nunits_i
1095 200 : xcoord(iunit) = r(1, iunit, imol)
1096 200 : ycoord(iunit) = r(2, iunit, imol)
1097 300 : zcoord(iunit) = r(3, iunit, imol)
1098 : END DO
1099 : EXIT
1100 : END IF
1101 : END DO
1102 300 : IF (lclus .EQV. .TRUE.) THEN
1103 600 : DO jmol = 1, nend
1104 500 : nunits_j = nunits(mol_type(jmol))
1105 600 : IF (clusmat(jmol) .EQ. 0 .AND. decision(jmol) .EQ. 0) THEN
1106 : !Calculating the distance between atoms
1107 600 : DO iunit = 1, nunits_i
1108 1400 : DO junit = 1, nunits_j
1109 800 : dx = xcoord(iunit) - r(1, junit, jmol)
1110 800 : dy = ycoord(iunit) - r(2, junit, jmol)
1111 800 : dz = zcoord(iunit) - r(3, junit, jmol)
1112 800 : dx = dx - abc(1)*ANINT(dx/abc(1))
1113 800 : dy = dy - abc(2)*ANINT(dy/abc(2))
1114 800 : dz = dz - abc(3)*ANINT(dz/abc(3))
1115 800 : rsquare = (dx*dx) + (dy*dy) + (dz*dz)
1116 : !Checking the distance based on rclus square(rclussq)
1117 1200 : IF (rsquare .LT. rclussquare) THEN
1118 0 : clusmat(jmol) = counter
1119 : END IF
1120 : END DO
1121 : END DO
1122 : END IF
1123 : END DO
1124 100 : DEALLOCATE (xcoord, ycoord, zcoord)
1125 : END IF
1126 : END DO
1127 : END DO
1128 :
1129 : !Putting cluster information in a cluster matrix
1130 20 : total_clus = counter
1131 :
1132 120 : DO imol = 1, counter
1133 620 : DO jmol = 1, nend
1134 600 : IF (imol .EQ. clusmat(jmol)) THEN
1135 100 : cluster(imol, jmol) = jmol
1136 : END IF
1137 : END DO
1138 : END DO
1139 20 : DEALLOCATE (r)
1140 20 : DEALLOCATE (decision)
1141 20 : DEALLOCATE (clusmat)
1142 :
1143 : ! end the timing
1144 20 : CALL timestop(handle)
1145 :
1146 60 : END SUBROUTINE cluster_search
1147 :
1148 : END MODULE mc_coordinates
1149 :
|
