Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 = 1.28558315_dp ! 0.6 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 32084 : DO imol = nstart, nend
137 24170 : IF (lall) jstart = imol + 1
138 24170 : nunits_i = nunits(mol_type(imol))
139 161714 : DO jmol = jstart, SUM(nchains(:))
140 93262 : IF (imol == jmol) CYCLE
141 91456 : nunits_j = nunits(mol_type(jmol))
142 :
143 353170 : DO iunit = 1, nunits_i
144 761138 : 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 430332 : (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 430332 : (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 430332 : (r(3, iunit, imol) - r(3, junit, jmol))/box_length(3))
155 :
156 430332 : dist = RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2
157 :
158 668612 : IF (dist < rmin) THEN
159 736 : loverlap = .TRUE.
160 736 : DEALLOCATE (r)
161 :
162 736 : CALL timestop(handle)
163 : RETURN
164 : END IF
165 :
166 : END DO
167 : END DO
168 : END DO
169 : END DO
170 :
171 7914 : DEALLOCATE (r)
172 :
173 : ! end the timing
174 7914 : 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 88 : ALLOCATE (boltz_weights(1:nswapmoves))
392 132 : ALLOCATE (loverlap_array(1:nswapmoves))
393 :
394 : ! initialize the arrays that need it
395 44 : loverlap_array(:) = .FALSE.
396 44 : loverlap = .FALSE.
397 44 : boltz_weights(:) = 0.0_dp
398 44 : trial_energy(:) = 0.0_dp
399 44 : 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 2 : boltz_weights(imove) = 0.0_dp
507 2 : CYCLE
508 : END IF
509 :
510 174 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
511 : CALL force_env_get(force_env, &
512 174 : potential_energy=bias_energy)
513 :
514 174 : trial_energy(imove) = (bias_energy - old_energy)
515 174 : exponent = -BETA*trial_energy(imove)
516 :
517 218 : IF (exponent > exp_max_val) THEN
518 0 : boltz_weights(imove) = max_val
519 174 : ELSE IF (exponent < exp_min_val) THEN
520 6 : 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 == 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 98 : all_overlaps = .FALSE.
556 98 : total_running_weight = total_running_weight + boltz_weights(imove)
557 98 : IF (total_running_weight >= rand*rosenbluth_weight) THEN
558 : choosen = imove
559 : EXIT
560 : END IF
561 : END DO
562 :
563 44 : 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) THEN
591 0 : CPABORT('CBMC swap move failed to select config')
592 : END IF
593 :
594 : ! if this is an old configuration, we always choose the first one
595 44 : IF (lremove) choosen = 1
596 :
597 : ! set the energy for the configuration
598 44 : choosen_energy = trial_energy(choosen)
599 :
600 : ! copy the coordinates to the force environment
601 1186 : DO iatom = 1, natoms_tot
602 4612 : particles%els(iatom)%r(1:3) = r(1:3, iatom, choosen)
603 : END DO
604 :
605 44 : DEALLOCATE (r_old)
606 44 : DEALLOCATE (r)
607 44 : DEALLOCATE (trial_energy)
608 44 : DEALLOCATE (boltz_weights)
609 44 : DEALLOCATE (loverlap_array)
610 :
611 : ! end the timing
612 44 : CALL timestop(handle)
613 :
614 44 : END SUBROUTINE generate_cbmc_swap_config
615 :
616 : ! **************************************************************************************************
617 : !> \brief rotates a molecule randomly around the center of mass,
618 : !> sequentially in x, y, and z directions
619 : !> \param r the coordinates of the molecule to rotate
620 : !> \param mass the mass of all the atoms in the molecule
621 : !> \param natoms the number of atoms in the molecule
622 : !> \param rng_stream the stream we pull random numbers from
623 : !>
624 : !> Use only in serial.
625 : !> \author MJM
626 : ! **************************************************************************************************
627 98 : SUBROUTINE rotate_molecule(r, mass, natoms, rng_stream)
628 :
629 : INTEGER, INTENT(IN) :: natoms
630 : REAL(KIND=dp), DIMENSION(1:natoms), INTENT(IN) :: mass
631 : REAL(KIND=dp), DIMENSION(1:3, 1:natoms), &
632 : INTENT(INOUT) :: r
633 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
634 :
635 : CHARACTER(len=*), PARAMETER :: routineN = 'rotate_molecule'
636 :
637 : INTEGER :: handle, iunit
638 : REAL(KIND=dp) :: cosdg, dgamma, rand, rx, rxnew, ry, &
639 : rynew, rz, rznew, sindg
640 : REAL(KIND=dp), DIMENSION(1:3) :: center_of_mass
641 :
642 : ! begin the timing of the subroutine
643 :
644 98 : CALL timeset(routineN, handle)
645 :
646 : ! find the center of mass of the molecule
647 98 : CALL get_center_of_mass(r(:, :), natoms, center_of_mass(:), mass(:))
648 :
649 : ! call a random number to figure out how far we're moving
650 98 : rand = rng_stream%next()
651 98 : dgamma = pi*(rand - 0.5E0_dp)*2.0E0_dp
652 :
653 : ! *** set up the rotation matrix ***
654 :
655 98 : cosdg = COS(dgamma)
656 98 : sindg = SIN(dgamma)
657 :
658 : ! *** ROTATE UNITS OF I AROUND X-AXIS ***
659 :
660 296 : DO iunit = 1, natoms
661 198 : ry = r(2, iunit) - center_of_mass(2)
662 198 : rz = r(3, iunit) - center_of_mass(3)
663 198 : rynew = cosdg*ry + sindg*rz
664 198 : rznew = cosdg*rz - sindg*ry
665 :
666 198 : r(2, iunit) = rynew + center_of_mass(2)
667 296 : r(3, iunit) = rznew + center_of_mass(3)
668 :
669 : END DO
670 :
671 : ! *** ROTATE UNITS OF I AROUND y-AXIS ***
672 :
673 296 : DO iunit = 1, natoms
674 198 : rx = r(1, iunit) - center_of_mass(1)
675 198 : rz = r(3, iunit) - center_of_mass(3)
676 198 : rxnew = cosdg*rx + sindg*rz
677 198 : rznew = cosdg*rz - sindg*rx
678 :
679 198 : r(1, iunit) = rxnew + center_of_mass(1)
680 296 : r(3, iunit) = rznew + center_of_mass(3)
681 :
682 : END DO
683 :
684 : ! *** ROTATE UNITS OF I AROUND z-AXIS ***
685 :
686 296 : DO iunit = 1, natoms
687 198 : rx = r(1, iunit) - center_of_mass(1)
688 198 : ry = r(2, iunit) - center_of_mass(2)
689 198 : rxnew = cosdg*rx + sindg*ry
690 198 : rynew = cosdg*ry - sindg*rx
691 :
692 198 : r(1, iunit) = rxnew + center_of_mass(1)
693 296 : r(2, iunit) = rynew + center_of_mass(2)
694 :
695 : END DO
696 :
697 : ! end the timing
698 98 : CALL timestop(handle)
699 :
700 98 : END SUBROUTINE rotate_molecule
701 :
702 : ! **************************************************************************************************
703 : !> \brief selects a molecule at random to perform a MC move on...you can specify
704 : !> the box the molecule should be in, its type, both, or neither
705 : !> \param mc_molecule_info the structure that contains some global molecule data
706 : !> \param start_atom the number of the first atom in the chosen molecule in relation
707 : !> to the force_env it's in
708 : !> \param box_number the box the chosen molecule is in
709 : !> \param molecule_type the type of molecule the chosen molecule is
710 : !> \param rng_stream the stream we pull random numbers from
711 : !> \param box if present, tells the routine which box to grab a molecule from
712 : !> \param molecule_type_old if present, tells the routine which molecule type to select from
713 : !> \author MJM
714 : ! **************************************************************************************************
715 816 : SUBROUTINE find_mc_test_molecule(mc_molecule_info, start_atom, &
716 : box_number, molecule_type, rng_stream, box, molecule_type_old)
717 :
718 : TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info
719 : INTEGER, INTENT(OUT) :: start_atom, box_number, molecule_type
720 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
721 : INTEGER, INTENT(IN), OPTIONAL :: box, molecule_type_old
722 :
723 : CHARACTER(LEN=*), PARAMETER :: routineN = 'find_mc_test_molecule'
724 :
725 : INTEGER :: handle, ibox, imol_type, imolecule, &
726 : jbox, molecule_number, nchains_tot, &
727 : start_mol
728 816 : INTEGER, DIMENSION(:), POINTER :: mol_type, nunits
729 816 : INTEGER, DIMENSION(:, :), POINTER :: nchains
730 : REAL(KIND=dp) :: rand
731 :
732 : ! begin the timing of the subroutine
733 :
734 816 : CALL timeset(routineN, handle)
735 :
736 816 : NULLIFY (nunits, mol_type, nchains)
737 : CALL get_mc_molecule_info(mc_molecule_info, nchains=nchains, nunits=nunits, &
738 816 : mol_type=mol_type)
739 :
740 : ! initialize the outgoing variables
741 816 : start_atom = 0
742 816 : box_number = 0
743 816 : molecule_type = 0
744 :
745 816 : IF (PRESENT(box) .AND. PRESENT(molecule_type_old)) THEN
746 : ! only need to find the atom number the molecule starts on
747 253 : rand = rng_stream%next()
748 253 : molecule_number = CEILING(rand*REAL(nchains(molecule_type_old, box), KIND=dp))
749 :
750 253 : start_mol = 1
751 253 : DO jbox = 1, box - 1
752 253 : start_mol = start_mol + SUM(nchains(:, jbox))
753 : END DO
754 :
755 : ! adjust to take into account molecules of other types in the box
756 253 : DO imol_type = 1, molecule_type_old - 1
757 253 : molecule_number = molecule_number + nchains(imol_type, box)
758 : END DO
759 :
760 253 : start_atom = 1
761 981 : DO imolecule = 1, molecule_number - 1
762 981 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
763 : END DO
764 :
765 563 : ELSE IF (PRESENT(box)) THEN
766 : ! any molecule in box...need to find molecule type and start atom
767 0 : rand = rng_stream%next()
768 0 : molecule_number = CEILING(rand*REAL(SUM(nchains(:, box)), KIND=dp))
769 :
770 0 : start_mol = 1
771 0 : DO jbox = 1, box - 1
772 0 : start_mol = start_mol + SUM(nchains(:, jbox))
773 : END DO
774 :
775 0 : molecule_type = mol_type(start_mol + molecule_number - 1)
776 :
777 : ! now the starting atom
778 0 : start_atom = 1
779 0 : DO imolecule = 1, molecule_number - 1
780 0 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
781 : END DO
782 :
783 563 : ELSE IF (PRESENT(molecule_type_old)) THEN
784 : ! any molecule of type molecule_type_old...need to find box number and start atom
785 563 : rand = rng_stream%next()
786 1162 : molecule_number = CEILING(rand*REAL(SUM(nchains(molecule_type_old, :)), KIND=dp))
787 :
788 : ! find which box it's in
789 563 : nchains_tot = 0
790 578 : DO ibox = 1, SIZE(nchains(molecule_type_old, :))
791 578 : IF (molecule_number <= nchains(molecule_type_old, ibox)) THEN
792 563 : box_number = ibox
793 563 : EXIT
794 : END IF
795 15 : molecule_number = molecule_number - nchains(molecule_type_old, ibox)
796 : END DO
797 :
798 563 : start_mol = 1
799 578 : DO jbox = 1, box_number - 1
800 608 : start_mol = start_mol + SUM(nchains(:, jbox))
801 : END DO
802 :
803 : ! now find the starting atom number
804 713 : DO imol_type = 1, molecule_type_old - 1
805 713 : molecule_number = molecule_number + nchains(imol_type, box_number)
806 : END DO
807 563 : start_atom = 1
808 3237 : DO imolecule = 1, molecule_number - 1
809 3237 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
810 : END DO
811 :
812 : ELSE
813 : ! no restrictions...need to find all pieces of data
814 0 : nchains_tot = 0
815 0 : DO ibox = 1, SIZE(nchains(1, :))
816 0 : nchains_tot = nchains_tot + SUM(nchains(:, ibox))
817 : END DO
818 0 : rand = rng_stream%next()
819 0 : molecule_number = CEILING(rand*REAL(nchains_tot, KIND=dp))
820 :
821 0 : molecule_type = mol_type(molecule_number)
822 :
823 : ! now which box it's in
824 0 : DO ibox = 1, SUM(nchains(1, :))
825 0 : IF (molecule_number <= SUM(nchains(:, ibox))) THEN
826 0 : box_number = ibox
827 0 : EXIT
828 : END IF
829 0 : molecule_number = molecule_number - SUM(nchains(:, ibox))
830 : END DO
831 :
832 : ! now find the starting atom number
833 0 : start_mol = 1
834 0 : DO jbox = 1, box_number - 1
835 0 : start_mol = start_mol + SUM(nchains(:, jbox))
836 : END DO
837 0 : start_atom = 1
838 0 : DO imolecule = 1, molecule_number - 1
839 0 : start_atom = start_atom + nunits(mol_type(start_mol + imolecule - 1))
840 : END DO
841 :
842 : END IF
843 :
844 : ! make sure things are good
845 816 : IF (PRESENT(box)) box_number = box
846 816 : IF (PRESENT(molecule_type_old)) molecule_type = molecule_type_old
847 :
848 816 : CPASSERT(start_atom /= 0)
849 816 : CPASSERT(box_number /= 0)
850 816 : CPASSERT(molecule_type /= 0)
851 :
852 : ! end the timing
853 816 : CALL timestop(handle)
854 :
855 816 : END SUBROUTINE find_mc_test_molecule
856 :
857 : ! **************************************************************************************************
858 : !> \brief generates an array that tells us which sides of the simulation
859 : !> cell we can increase or decrease using a discrete volume move
860 : !> \param cell the lengths of the sides of the cell
861 : !> \param discrete_array the array that indicates which sides we can move
862 : !> \param step_size the size of the discrete volume move
863 : !>
864 : !> Suitable for parallel.
865 : !> \author MJM
866 : ! **************************************************************************************************
867 0 : SUBROUTINE create_discrete_array(cell, discrete_array, step_size)
868 :
869 : ! 1 is for increase, 2 is for decrease
870 : ! 1 is for "yes, we can do the move", 0 is for no
871 :
872 : REAL(dp), DIMENSION(1:3), INTENT(IN) :: cell
873 : INTEGER, DIMENSION(1:3, 1:2), INTENT(OUT) :: discrete_array
874 : REAL(dp), INTENT(IN) :: step_size
875 :
876 : INTEGER :: iside
877 : REAL(dp) :: high_value, length1, length2, low_value
878 :
879 0 : discrete_array(:, :) = 0
880 :
881 0 : length1 = ABS(cell(1) - cell(2))
882 0 : length2 = ABS(cell(2) - cell(3))
883 :
884 : ! now let's figure out all the different cases
885 0 : IF (length1 < 0.01_dp*step_size .AND. &
886 : length2 < 0.01_dp*step_size) THEN
887 : ! all sides are equal, so we can move up or down
888 0 : discrete_array(1:3, 1) = 1
889 0 : discrete_array(1:3, 2) = 1
890 : ELSE
891 :
892 : ! find the low value and the high value
893 0 : high_value = -1.0_dp
894 0 : low_value = cell(1)*cell(2)*cell(3)
895 0 : DO iside = 1, 3
896 0 : IF (cell(iside) < low_value) low_value = cell(iside)
897 0 : IF (cell(iside) > high_value) high_value = cell(iside)
898 : END DO
899 0 : DO iside = 1, 3
900 : ! now we see if the value is a high value or a low value...it can only be
901 : ! one of the two
902 0 : IF (ABS(cell(iside) - low_value) < 0.01_dp*step_size) THEN
903 : ! low value, we can only increase the cell size
904 0 : discrete_array(iside, 1) = 1
905 0 : discrete_array(iside, 2) = 0
906 : ELSE
907 : ! high value, we can only decrease the cell size
908 0 : discrete_array(iside, 1) = 0
909 0 : discrete_array(iside, 2) = 1
910 : END IF
911 : END DO
912 : END IF
913 :
914 0 : END SUBROUTINE create_discrete_array
915 :
916 : ! **************************************************************************************************
917 : !> \brief generates an insertion point in either the "in" or the "out" volume
918 : !> of a target atom, where the "in" volume is a shell with inner radius
919 : !> rmin and outer radius rmax
920 : !> \param rmin the minimum AVBMC radius for the shell around the target
921 : !> \param rmax the maximum AVBMC radius for the shell around the target
922 : !> \param r_target the coordinates of the target atom
923 : !> \param move_type generate configs in the "in" or "out" volume
924 : !> \param r_insert the output insertion site
925 : !> \param abc the lengths of the sides of the simulation box
926 : !> \param rng_stream the random number stream that we draw from
927 : !>
928 : !> Use only in serial.
929 : !> \author MJM
930 : ! **************************************************************************************************
931 0 : SUBROUTINE generate_avbmc_insertion(rmin, rmax, r_target, &
932 : move_type, r_insert, abc, rng_stream)
933 :
934 : REAL(KIND=dp), INTENT(IN) :: rmin, rmax
935 : REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: r_target
936 : CHARACTER(LEN=*), INTENT(IN) :: move_type
937 : REAL(KIND=dp), DIMENSION(1:3), INTENT(OUT) :: r_insert
938 : REAL(KIND=dp), DIMENSION(1:3), INTENT(IN) :: abc
939 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
940 :
941 : INTEGER :: i
942 : REAL(dp) :: dist, eta_1, eta_2, eta_sq, rand
943 : REAL(dp), DIMENSION(1:3) :: RIJ
944 :
945 0 : r_insert(1:3) = 0.0_dp
946 :
947 0 : IF (move_type == 'in') THEN
948 : ! generate a random unit vector, from Allen and Tildesley
949 : DO
950 0 : eta_1 = rng_stream%next()
951 0 : eta_2 = rng_stream%next()
952 0 : eta_sq = eta_1**2 + eta_2**2
953 0 : IF (eta_sq < 1.0_dp) THEN
954 0 : r_insert(1) = 2.0_dp*eta_1*SQRT(1.0_dp - eta_sq)
955 0 : r_insert(2) = 2.0_dp*eta_2*SQRT(1.0_dp - eta_sq)
956 0 : r_insert(3) = 1.0_dp - 2.0_dp*eta_sq
957 : EXIT
958 : END IF
959 : END DO
960 :
961 : ! now scale that vector to be within the "in" region
962 0 : rand = rng_stream%next()
963 : r_insert(1:3) = r_insert(1:3)*(rand*(rmax**3 - rmin**3) + rmin**3)** &
964 0 : (1.0_dp/3.0_dp)
965 :
966 0 : r_insert(1:3) = r_target(1:3) + r_insert(1:3)
967 : ELSE
968 :
969 : ! find a new insertion point somewhere in the box
970 : DO
971 0 : DO i = 1, 3
972 0 : rand = rng_stream%next()
973 0 : r_insert(i) = rand*abc(i)
974 : END DO
975 :
976 : ! make sure it's not in the "in" region
977 : RIJ(1) = r_insert(1) - r_target(1) - abc(1)* &
978 0 : ANINT((r_insert(1) - r_target(1))/abc(1))
979 : RIJ(2) = r_insert(2) - r_target(2) - abc(2)* &
980 0 : ANINT((r_insert(2) - r_target(2))/abc(2))
981 : RIJ(3) = r_insert(3) - r_target(3) - abc(3)* &
982 0 : ANINT((r_insert(3) - r_target(3))/abc(3))
983 :
984 0 : dist = RIJ(1)**2 + RIJ(2)**2 + RIJ(3)**2
985 :
986 0 : IF (dist < rmin**2 .OR. dist > rmax**2) THEN
987 : EXIT
988 : END IF
989 :
990 : END DO
991 : END IF
992 :
993 0 : END SUBROUTINE generate_avbmc_insertion
994 :
995 : ! *****************************************************************************
996 : ! **************************************************************************************************
997 : !> \brief determine the number of cluster present in the given configuration
998 : !> based on the rclus value
999 : !> \param mc_par the mc parameters for the force env
1000 : !> \param force_env the force environment containing the coordinates
1001 : !> \param cluster ...
1002 : !> \param nchains the number of molecules of each type in the box
1003 : !> \param nunits the number of interaction sites for each molecule
1004 : !> \param mol_type an array that indicates the type of each molecule
1005 : !> \param total_clus ...
1006 : !> \par
1007 : !> Original Multiparticle/Cluster Translation paper:
1008 : !> Orkoulas, Gerassimos, and Athanassios Z. Panagiotopoulos. Free energy and
1009 : !> phase equilibria for the restricted primitive model of ionic fluids from Monte
1010 : !> Carlo simulations. J. Chem. Phys. 1994,101.2,1452-1459.
1011 : !> \author Himanshu Goel
1012 : ! **************************************************************************************************
1013 :
1014 20 : SUBROUTINE cluster_search(mc_par, force_env, cluster, nchains, nunits, mol_type, total_clus)
1015 :
1016 : TYPE(mc_simpar_type), POINTER :: mc_par
1017 : TYPE(force_env_type), POINTER :: force_env
1018 : INTEGER, DIMENSION(:, :), INTENT(INOUT) :: cluster
1019 : INTEGER, DIMENSION(:), INTENT(IN) :: nchains, nunits, mol_type
1020 : INTEGER, INTENT(INOUT) :: total_clus
1021 :
1022 : CHARACTER(len=*), PARAMETER :: routineN = 'cluster_search'
1023 :
1024 : INTEGER :: counter, handle, imol, iunit, jmol, &
1025 : junit, nend, nstart, nunit, nunits_i, &
1026 : nunits_j
1027 20 : INTEGER, ALLOCATABLE, DIMENSION(:) :: clusmat, decision
1028 : LOGICAL :: lclus
1029 : REAL(KIND=dp) :: dx, dy, dz, rclus, rclussquare, rsquare
1030 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xcoord, ycoord, zcoord
1031 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: r
1032 : REAL(KIND=dp), DIMENSION(1:3) :: abc
1033 : TYPE(cell_type), POINTER :: cell
1034 : TYPE(cp_subsys_type), POINTER :: oldsys
1035 : TYPE(particle_list_type), POINTER :: particles
1036 :
1037 : ! begin the timing of the subroutine
1038 :
1039 20 : CALL timeset(routineN, handle)
1040 :
1041 20 : NULLIFY (oldsys, particles)
1042 :
1043 : ! Getting Particle Coordinates
1044 :
1045 20 : CALL force_env_get(force_env, cell=cell, subsys=oldsys)
1046 20 : CALL get_cell(cell, abc=abc)
1047 20 : CALL cp_subsys_get(oldsys, particles=particles)
1048 20 : CALL get_mc_par(mc_par, rclus=rclus)
1049 :
1050 120 : ALLOCATE (r(1:3, 1:MAXVAL(nunits), 1:SUM(nchains)))
1051 :
1052 : ! Arranging particles coordinates into an easier matrix
1053 40 : nend = SUM(nchains(:))
1054 : junit = 0
1055 120 : DO imol = 1, nend
1056 100 : nunit = nunits(mol_type(imol))
1057 320 : DO iunit = 1, nunit
1058 200 : junit = junit + 1
1059 900 : r(1:3, iunit, imol) = particles%els(junit)%r(1:3)
1060 : END DO
1061 : END DO
1062 :
1063 20 : counter = 0
1064 :
1065 : ! Allocating the size of matrix and decision matrix
1066 80 : ALLOCATE (clusmat(nend), decision(nend))
1067 :
1068 : !Initialize
1069 120 : DO imol = 1, nend
1070 100 : decision(imol) = 0
1071 120 : clusmat(imol) = 0
1072 : END DO
1073 :
1074 20 : rclussquare = rclus*rclus
1075 : ! Starting the cluster count loop
1076 720 : DO WHILE (SUM(decision) < nend)
1077 300 : DO nstart = 1, nend
1078 300 : IF (clusmat(nstart) == 0) THEN
1079 100 : counter = counter + 1
1080 100 : clusmat(nstart) = counter
1081 100 : EXIT
1082 : END IF
1083 : END DO
1084 :
1085 100 : lclus = .TRUE.
1086 220 : DO WHILE (lclus .EQV. .TRUE.)
1087 900 : DO imol = 1, nend
1088 800 : nunits_i = nunits(mol_type(imol))
1089 : ! Allocating the xcoord,ycoord,zcoord based upon the size of molecule nunits
1090 800 : lclus = .FALSE.
1091 900 : IF (clusmat(imol) == counter .AND. decision(imol) == 0) THEN
1092 500 : ALLOCATE (xcoord(nunits_i), ycoord(nunits_i), zcoord(nunits_i))
1093 100 : decision(imol) = 1
1094 100 : lclus = .TRUE.
1095 300 : DO iunit = 1, nunits_i
1096 200 : xcoord(iunit) = r(1, iunit, imol)
1097 200 : ycoord(iunit) = r(2, iunit, imol)
1098 300 : zcoord(iunit) = r(3, iunit, imol)
1099 : END DO
1100 : EXIT
1101 : END IF
1102 : END DO
1103 200 : IF (lclus .EQV. .TRUE.) THEN
1104 600 : DO jmol = 1, nend
1105 500 : nunits_j = nunits(mol_type(jmol))
1106 600 : IF (clusmat(jmol) == 0 .AND. decision(jmol) == 0) THEN
1107 : !Calculating the distance between atoms
1108 600 : DO iunit = 1, nunits_i
1109 1400 : DO junit = 1, nunits_j
1110 800 : dx = xcoord(iunit) - r(1, junit, jmol)
1111 800 : dy = ycoord(iunit) - r(2, junit, jmol)
1112 800 : dz = zcoord(iunit) - r(3, junit, jmol)
1113 800 : dx = dx - abc(1)*ANINT(dx/abc(1))
1114 800 : dy = dy - abc(2)*ANINT(dy/abc(2))
1115 800 : dz = dz - abc(3)*ANINT(dz/abc(3))
1116 800 : rsquare = (dx*dx) + (dy*dy) + (dz*dz)
1117 : !Checking the distance based on rclus square(rclussq)
1118 1200 : IF (rsquare < rclussquare) THEN
1119 0 : clusmat(jmol) = counter
1120 : END IF
1121 : END DO
1122 : END DO
1123 : END IF
1124 : END DO
1125 100 : DEALLOCATE (xcoord, ycoord, zcoord)
1126 : END IF
1127 : END DO
1128 : END DO
1129 :
1130 : !Putting cluster information in a cluster matrix
1131 20 : total_clus = counter
1132 :
1133 120 : DO imol = 1, counter
1134 620 : DO jmol = 1, nend
1135 600 : IF (imol == clusmat(jmol)) THEN
1136 100 : cluster(imol, jmol) = jmol
1137 : END IF
1138 : END DO
1139 : END DO
1140 20 : DEALLOCATE (r)
1141 20 : DEALLOCATE (decision)
1142 20 : DEALLOCATE (clusmat)
1143 :
1144 : ! end the timing
1145 20 : CALL timestop(handle)
1146 :
1147 60 : END SUBROUTINE cluster_search
1148 :
1149 : END MODULE mc_coordinates
1150 :
|
