chemshell_qchem / blob
1 ! **********************************************************************
2 ! ** **
3 ! ** Parallel optimisation **
4 ! ** **
5 ! **********************************************************************
6 !!****h* DL-FIND/parallel_opt
7 !!
8 !! NAME
9 !! parallel_opt
10 !!
11 !! FUNCTION
12 !! Performs a parallel optimisation; currently either stochastic search
13 !! or genetic algorithm.
14 !!
15 !! COPYRIGHT
16 !!
17 !! Copyright 2007 Johannes Kaestner (kaestner@theochem.uni-stuttgart.de),
18 !! Tom Keal (thomas.keal@stfc.ac.uk),
19 !! Joanne Carr (j.m.carr@dl.ac.uk)
20 !!
21 !! This file is part of DL-FIND.
22 !!
23 !! DL-FIND is free software: you can redistribute it and/or modify
24 !! it under the terms of the GNU Lesser General Public License as
25 !! published by the Free Software Foundation, either version 3 of the
26 !! License, or (at your option) any later version.
27 !!
28 !! DL-FIND is distributed in the hope that it will be useful,
29 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
30 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 !! GNU Lesser General Public License for more details.
32 !!
33 !! You should have received a copy of the GNU Lesser General Public
34 !! License along with DL-FIND. If not, see
35 !! <http://www.gnu.org/licenses/>.
36 !!
37 !! SOURCE
38 !!****
40 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41 !!****f* parallel_opt/dlf_parallel_opt
42 !!
43 !! FUNCTION
44 !! Set up, perform and close down after, the parallel optimisation
45 !! Note that two numbers are hardwired in the code at present --
46 !! small, mainly used as the tolerance on fractional energy differences to
47 !! distinguish between structures; and the reduction of the radius by a
48 !! factor of 4 after the initial population is set up in the genetic algorithm
49 !! (search for the comment: Scale the radius from its initial value).
50 !!
51 !! INPUTS
52 !!
53 !! glob%xcoords
54 !!
55 !! for SS --
56 !! glob%po_pop_size
57 !! glob%po_radius(:)
58 !! glob%po_contraction
59 !! glob%po_tolerance_r(:)
60 !! glob%po_tolerance_g
61 !! glob%po_distribution
62 !! glob%po_maxcycle
63 !! glob%po_scalefac
64 !!
65 !! for GA --
66 !! glob%po_pop_size
67 !! glob%po_radius(:)
68 !! glob%po_tolerance_g
69 !! glob%po_maxcycle
70 !! glob%po_init_pop_size
71 !! glob%po_reset
72 !! glob%po_mutation_rate
73 !! glob%po_death_rate
74 !! glob%po_nsave
75 !!
76 !! OUTPUTS
77 !!
78 !! tconv
79 !!
80 !! Local variables:
81 !! * xcoords_best(:,:)
82 !! * energy_best
83 !! * for GA, if nsave/=0 -- energies_save(:)
84 !! nevals_save(:)
85 !! xcoords_save(:,:,:)
86 !!
87 !! Files, if converged -- best.xyz
88 !! best_active.xyz
89 !!
90 !! SYNOPSIS
91 subroutine dlf_parallel_opt(trestarted_report, tconv &
92 #ifdef GAMESS
93 ,core&
94 #endif
95 )
96 !! SOURCE
97 use dlf_parameter_module, only: rk
98 use dlf_global, only: glob,stderr,stdout,printl,printf
99 use dlf_stat, only: stat
100 use dlf_allocate, only: allocate,deallocate
101 use dlf_checkpoint
102 use dlf_sort_module
103 implicit none
104 #ifdef GAMESS
105 real(rk) :: core(*) ! GAMESS memory, not used in DL-FIND
106 #endif
107 integer :: status, iimage, i, j, k, seed_size, nmutations
108 integer :: nresets, noffspring, lower_index, arr_size
109 integer,dimension(1) :: position
110 integer,dimension(8) :: time_data
111 integer,allocatable :: seed(:), nevals_save(:)
112 logical :: stochastic, genetic, trestarted_report, tconv, texit
113 logical :: dummy_logical, trestarted, tno_diversity
114 real(rk) :: random, energy_best
115 real(rk) :: mean_init_e, mean_init_e2, sigma_init_e
116 real(rk),allocatable :: pop_energies(:), init_pop_energies(:)
117 real(rk),allocatable :: icoords_best(:), igradient_best(:), xcoords_best(:,:)
118 real(rk),allocatable :: energies_save(:), xcoords_save(:,:,:)
119 real(rk),allocatable :: pop_icoords(:,:), init_pop_icoords(:,:)
120 real(rk),allocatable :: pop_xgradient(:,:,:)
122 ! may want to change this for different systems...
123 real(rk), parameter :: small=1.0D-12
125 ! **********************************************************************
127 ! Initialise variables
128 stochastic = (glob%iopt == 51)
129 genetic = (glob%iopt == 52)
130 tconv = .false.
131 texit = .false.
132 seed_size = 1
133 energy_best = huge(1.0D0)
134 iimage = 1
135 ! stat%sene = 0
136 ! stat%ccycle = 0
137 ! GA variables
138 if (genetic) then
139 nresets = 0
140 tno_diversity = .false.
142 if (glob%po_pop_size < 4) then
143 write(stdout,'(1x,a)')"ERROR: working population size is too small"
144 write(stdout,'(1x,a)')"Size must be >= 4 for genetic algorithm"
145 call dlf_fail("Input population size too small")
146 end if
148 if (glob%po_init_pop_size < glob%po_pop_size) then
149 write(stdout,'(1x,a)')"ERROR: initial population size is too small"
150 write(stdout,'(1x,a)')"Size must be >= working population size"
151 call dlf_fail("Input initial population size too small")
152 end if
154 nmutations = nint(glob%po_pop_size * glob%nivar * glob%po_mutation_rate)
155 ! use nint above to allow nmutations = 0
156 if (nmutations > glob%po_pop_size * glob%nivar) &
157 &nmutations = glob%po_pop_size * glob%nivar
158 noffspring = int(glob%po_pop_size * glob%po_death_rate) + 1
159 ! int truncates, so noffspring will never be zero
160 ! check noffspring against the upper limit of glob%po_pop_size - 2
161 ! (want to leave at least two individuals as possible parents)
162 if (noffspring > glob%po_pop_size - 2) noffspring = glob%po_pop_size - 2
163 if (mod(noffspring,2) == 1) then
164 ! want to replace the two parents with their two offspring, so ensure
165 ! noffspring is even by increasing the value by 1 if it is odd
166 noffspring = noffspring + 1
167 if (noffspring > glob%po_pop_size - 2) noffspring = noffspring - 2
168 end if
170 write(stdout,'(1x,a,i10)')"In dlf_parallel_opt, noffspring = ",noffspring
171 write(stdout,'(1x,a,i10)')"In dlf_parallel_opt, nmutations = ",nmutations
172 end if
174 ! Initialise random number generator on the zeroth processor only
175 ! (no other processor should generate any random numbers...)
176 ! (pgf95 is rather picky about the random number seed, hence the
177 ! complicated procedure below...)
179 if (glob%iam == 0) then
180 call random_seed(size = seed_size)
181 call allocate(seed, seed_size)
182 seed = 1
183 call date_and_time(values = time_data)
184 seed(1) = time_data(8)
185 do i = 1, seed_size - 1
186 seed(i + 1) = time_data(8 - mod(i,8)) * seed(1) * i
187 ! use smallest units of time_data first, and make the repeats different
188 if (seed(i + 1) == 0) seed(i + 1) = i
189 end do
190 seed = 11*seed ! so seed does not just contain small values
191 ! Also helps to ensure that not all the elements are even
192 ! (which hinders performance for pgf95)
193 call random_seed(put = seed)
194 !!! for testing call random_seed(get = seed)
195 write(stdout,'(1x,a)')"Random number seed = "
196 do i = 1, seed_size
197 write(stdout,*) seed(i)
198 end do
199 write(stdout,'(1x,a)')"End of random number seed"
200 call random_number(random)
201 write(stdout,*)"random 1 ",random
202 call random_number(random)
203 write(stdout,*)"random 2 ",random
205 call deallocate(seed)
206 end if
208 ! allocate storage
209 call allocate(pop_icoords, glob%po_pop_size, glob%nivar)
210 call allocate(pop_xgradient, glob%po_pop_size, 3, glob%nat)
211 call allocate(pop_energies, glob%po_pop_size)
212 call allocate(icoords_best, glob%nivar)
213 call allocate(igradient_best, glob%nivar)
214 call allocate(xcoords_best, 3, glob%nat)
215 call allocate(glob%po_radius, glob%nivar)
216 call allocate(glob%po_tolerance_r, glob%nivar)
218 ! Set the sample and tolerance radii arrays
219 ! Note that on restart, the base values glob%po_radius_base and glob%po_tol_r_base
220 ! are obtained from the restart file, but the contents of glob%po_radius_scaling(:)
221 ! must still be determined from the passed tmpcoords2 array in dlf_get_params.
222 ! glob%po_radius_scaling(:) can therefore be different from the original if desired.
223 arr_size = SIZE(glob%po_radius_scaling, 1)
224 if (arr_size == 1 .and. glob%nivar > 1) then
225 glob%po_radius(:) = glob%po_radius_base
226 glob%po_tolerance_r(:) = glob%po_tol_r_base
227 elseif (arr_size == glob%nivar) then
228 glob%po_radius(:) = glob%po_radius_base * glob%po_radius_scaling(:)
229 glob%po_tolerance_r(:) = glob%po_tol_r_base * glob%po_radius_scaling(:)
230 else
231 call dlf_fail("Radii scaling factors incompatible with working coordinates")
232 end if
234 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
235 ! Setup of initial population or individual
236 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238 if (stochastic) then
240 lower_index = 1
242 ! read checkpoint file OR sort out initial individual
243 if (glob%restart == 1) then
245 call clock_start("CHECKPOINT")
246 call dlf_checkpoint_po_read
247 call clock_stop("CHECKPOINT")
248 if (.not.trestarted) call dlf_fail("Restart attempt failed")
250 else
252 ! evaluate initial point
253 call clock_start("EANDG")
254 call dlf_get_gradient(glob%nvar,glob%xcoords,energy_best, &
255 glob%xgradient,iimage,&
256 #ifdef GAMESS
257 core,&
258 #endif
259 status)
260 stat%sene = stat%sene + 1
261 stat%pene = stat%pene + 1
262 call clock_stop("EANDG")
264 ! check of NaN in the energy (comparison of NaN with any number
265 ! returns .false. , pgf90 does not understand isnan() )
266 if ( abs(energy_best) > huge(1.D0) ) then
267 status = 1
268 else
269 if (.not. abs(energy_best) < huge(1.D0) ) status = 1
270 end if
272 if (status/=0) then
273 call dlf_report(trestarted_report)
274 call dlf_fail("Energy evaluation failed")
275 end if
277 xcoords_best(:,:) = glob%xcoords(:,:)
279 ! convert glob%xcoords and glob%xgradient from Cartesians to the
280 ! working (i-labelled) coordinates glob%icoords and glob%igradient
281 call clock_start("COORDS")
282 call dlf_coords_xtoi(dummy_logical,dummy_logical,iimage,0)
283 call clock_stop("COORDS")
285 icoords_best(:) = glob%icoords(:)
286 igradient_best(:) = glob%igradient(:)
288 end if
290 else if (genetic) then
292 lower_index = 2
293 call allocate(init_pop_icoords,glob%po_init_pop_size,glob%nivar)
294 call allocate(init_pop_energies,glob%po_init_pop_size)
296 if (glob%po_nsave > 0) then
297 call allocate(xcoords_save,glob%po_nsave,3,glob%nat)
298 call allocate(energies_save,glob%po_nsave)
299 call allocate(nevals_save,glob%po_nsave)
300 energies_save(:) = huge(1.0D0)
301 nevals_save(:) = 0
302 xcoords_save(:,:,:) = 0.0D0
303 end if
305 ! read checkpoint file OR sort out initial population
306 if (glob%restart == 1) then
308 call clock_start("CHECKPOINT")
309 call dlf_checkpoint_po_read
310 call clock_stop("CHECKPOINT")
311 if (.not.trestarted) call dlf_fail("Restart attempt failed")
313 else
315 ! convert glob%xcoords from Cartesians to the
316 ! working (i-labelled) coordinates glob%icoords
317 ! The glob%[x,i]gradient arrays contain junk at this point
318 glob%xgradient = 0.0D0 ! Safest to initialize though
319 call clock_start("COORDS")
320 call dlf_coords_xtoi(dummy_logical,dummy_logical,iimage,0)
321 call clock_stop("COORDS")
323 ! Make initial population...
324 if (glob%iam == 0) then
325 call clock_start("FORMSTEP")
326 call dlf_genetic_initialpop(init_pop_icoords(:,:),&
327 & glob%po_init_pop_size, glob%icoords(:))
328 call clock_stop("FORMSTEP")
329 end if
331 call dlf_global_real_bcast(init_pop_icoords(:,:), &
332 & glob%po_init_pop_size*glob%nivar, 0)
334 ! ...then work out the energy of each individual...
335 init_pop_energies(:) = 0.0D0
336 do j = 1, glob%po_init_pop_size
337 k = mod(j,glob%ntasks)
338 stat%sene = stat%sene + 1
339 if (k == glob%mytask) then
341 glob%icoords(:) = init_pop_icoords(j,:)
343 ! convert glob%icoords (used as a temporary storage array) back
344 ! to Cartesians in glob%xcoords, before call to dlf_get_gradient
346 call clock_start("COORDS")
347 call dlf_coords_itox(iimage)
348 call clock_stop("COORDS")
350 call clock_start("EANDG")
351 !!! using pop_xgradient array as a dummy in the call, because IN THE CURRENT
352 !!! CODE it will not be used before it is filled with the proper values, and
353 !!! to save memory I don't want an unnecessary init_pop_xgradient array.
354 call dlf_get_gradient(glob%nvar,glob%xcoords, &
355 init_pop_energies(j), pop_xgradient(1,:,:),iimage,&
356 #ifdef GAMESS
357 core,&
358 #endif
359 status)
360 stat%pene = stat%pene + 1
361 call clock_stop("EANDG")
363 ! check of NaN in the energy (comparison of NaN with any number
364 ! returns .false. , pgf90 does not understand isnan() )
365 if ( abs(init_pop_energies(j)) > huge(1.D0) ) then
366 status = 1
367 else
368 if (.not. abs(init_pop_energies(j)) < huge(1.D0) ) status = 1
369 end if
371 if (status/=0) then
372 init_pop_energies(j) = huge(1.0D0)
373 !pop_xgradient(j,:,:) = huge(1.0D0)
374 end if
376 endif
377 enddo
379 !!!call dlf_global_real_sum(init_pop_energies(:), glob%po_init_pop_size)
380 !!!init_pop_energies(:) = init_pop_energies(:)/glob%nprocs_per_task
381 call dlf_tasks_real_sum(init_pop_energies(:), glob%po_init_pop_size)
383 ! ...then sort the population on energy...
384 call dlf_sort(init_pop_icoords(:,:), init_pop_energies(:), &
385 & glob%po_init_pop_size, glob%nivar)
387 ! Write information for the starting population
388 mean_init_e = sum(init_pop_energies(:)) / dble(glob%po_init_pop_size)
389 mean_init_e2 =sum(init_pop_energies(:)**2)/dble(glob%po_init_pop_size)
390 sigma_init_e = sqrt(mean_init_e2 - mean_init_e**2)
392 write(stdout,'(1x,a)')"Initial population:"
393 write(stdout,'(1x,a,es16.9)')"Lowest energy = ",init_pop_energies(1)
394 write(stdout,'(1x,a,es16.9)')"Mean energy = ",mean_init_e
395 write(stdout,'(1x,a,es16.9)')"Standard deviation of the energy = ",&
396 &sigma_init_e
398 ! ...and put the glob%po_pop_size fittest (=lowest-energy) individuals
399 ! into the working population, pop_icoords(:,:) and pop_energies(:)
400 pop_icoords(:,:) = init_pop_icoords(:glob%po_pop_size,:)
401 pop_energies(:) = init_pop_energies(:glob%po_pop_size)
403 end if
405 ! Scale the radius from its initial value, for use in setting up the
406 ! population, to something more appropriate for the mutations.
408 glob%po_radius(:) = 0.25D0*glob%po_radius(:)
410 else
412 call dlf_fail("Parallel optimisation option requested that does not &
413 &yet exist")
415 end if
417 call flush(stdout)
419 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420 ! End setup of initial population or individual
421 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
424 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
425 ! Main cycles loop
426 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
428 do i = 1, glob%po_maxcycle
430 stat%ccycle = stat%ccycle + 1
432 ! Generate the individuals in the population for this cycle and store them in
433 ! pop_icoords(:,:)
434 if (glob%iam == 0) then
435 call clock_start("FORMSTEP")
436 if (stochastic) then
437 call dlf_stoch_makepop
438 else if (genetic) then
439 call dlf_genetic_makepop
440 end if
441 call clock_stop("FORMSTEP")
442 end if
444 if (genetic) then
445 call dlf_global_log_bcast(tno_diversity, 1, 0)
446 if (tno_diversity) then
447 write(stdout,'(1x,a)')"dlf_parallel_opt WARNING"
448 write(stdout,'(1x,a)')"Insufficient genetic diversity in population"
449 !exit ! don't necesarily need to exit here
450 write(stdout,'(1x,a)')"Resetting population from known individual &
451 &of lowest energy"
452 write(stdout,'(1x,a)')"(Could also try a higher mutation rate)"
453 if (glob%iam == 0) then
454 call clock_start("FORMSTEP")
455 call dlf_genetic_initialpop(pop_icoords(:,:), &
456 & glob%po_pop_size, pop_icoords(1,:))
457 ! Note, dlf_genetic_initialpop doesn't calculate energies
458 call clock_stop("FORMSTEP")
459 end if
460 tno_diversity = .false.
461 else
462 if (nmutations > 0) then
463 if (glob%iam == 0) call dlf_genetic_mutate
464 end if
465 end if
466 end if
468 ! This call can go here rather than before the if (genetic) block above
469 ! because the subroutine calls that require pop_icoords(:,:) to be known
470 ! should only be made from the master processor (glob%iam ==0), which knows
471 ! the full pop_icoords after the call to dlf_genetic_makepop
473 call dlf_global_real_bcast(pop_icoords(:,:),glob%po_pop_size*glob%nivar,0)
475 ! This (re)initialization is necessary for the parallel implementation,
476 ! when MPI_Allreduce calls are made for the pop_energies(:) and
477 ! pop_xgradient(:,:,:) arrays.
478 call dlf_init_engarrays
480 ! calculate energy and gradient for each individual in the population
481 ! (apart from the first individual in a GA run...)
482 call dlf_get_engarrays
484 ! Make the full pop_energies (and pop_xgradient) arrays known on all processors
485 !!!call dlf_global_real_sum(pop_energies(:), glob%po_pop_size)
486 !!!pop_energies(:) = pop_energies(:)/glob%nprocs_per_task
487 !!!call dlf_global_real_sum(pop_xgradient(:,:,:), glob%po_pop_size*glob%nvar)
488 !!!pop_xgradient(:,:,:) = pop_xgradient(:,:,:)/glob%nprocs_per_task
489 call dlf_tasks_real_sum(pop_energies(:), glob%po_pop_size)
490 call dlf_tasks_real_sum(pop_xgradient(:,:,:), glob%po_pop_size*glob%nvar)
492 ! Prepare for next cycle; test for convergence or failure
493 if (stochastic) then
494 call dlf_stoch_endcycle
495 if (texit) exit
496 else if (genetic) then
497 call dlf_sort(pop_icoords(:,:), pop_xgradient(:,:,:), &
498 & pop_energies(:), glob%po_pop_size, glob%nivar, 3, glob%nat)
499 call dlf_genetic_endcycle
500 end if
502 if ( (mod(i,glob%po_reset) == 0) .and. genetic) then
503 nresets = nresets + 1
504 write(stdout,'(1x,a)')"Resetting population from known individual of &
505 &lowest energy"
507 if (glob%iam == 0) then
508 call clock_start("FORMSTEP")
509 call dlf_genetic_initialpop(pop_icoords(:,:), &
510 & glob%po_pop_size, pop_icoords(1,:))
511 call clock_stop("FORMSTEP")
512 end if
514 call dlf_global_real_bcast(pop_icoords(:,:),&
515 &glob%po_pop_size*glob%nivar,0)
516 call dlf_init_engarrays
517 call dlf_get_engarrays
518 call dlf_tasks_real_sum(pop_energies(:), glob%po_pop_size)
519 call dlf_tasks_real_sum(pop_xgradient(:,:,:), &
520 &glob%po_pop_size*glob%nvar)
521 call dlf_sort(pop_icoords(:,:), pop_xgradient(:,:,:), &
522 & pop_energies(:), glob%po_pop_size, glob%nivar, 3, glob%nat)
524 end if
526 call flush(stdout)
528 end do ! the cycles loop
530 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
531 ! End of main cycles loop
532 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
534 if (genetic .and. tconv .and. glob%po_nsave > 0) then
535 write(stdout,'(1x,a)')"End of genetic algorithm"
536 write(stdout,'(1x,a,i10,a)')"List of ",glob%po_nsave," lowest-energy minima"
537 do i = 1, glob%po_nsave
538 if (energies_save(i) == huge(1.0D0)) then
539 write(stdout,'(1x,a)')"Warning: no more unique minima found"
540 exit
541 end if
542 write(stdout,'(1x,i10)') i
543 write(stdout,'(1x,a,es16.9)')"Energy = ",energies_save(i)
544 write(stdout,'(1x,a,i10,a)')"Found after ",nevals_save(i),&
545 &" energy evaluations"
546 write(stdout,'(1x,a)')"Cartesian coordinates:"
547 write(stdout,'(3(1x,es16.9))') xcoords_save(i,:,:)
548 call dlf_put_coords(glob%nvar,-i,energies_save(i),xcoords_save(i,:,:)&
549 &,glob%iam)
550 end do
551 end if
553 call flush(stdout)
554 call dlf_po_destroy ! deallocate memory
555 return
557 contains
558 !!****
560 ! Internal procedures as they should not be called from any subroutine apart
561 ! from dlf_parallel_opt
564 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
565 !!****f* parallel_opt/dlf_init_engarrays
566 !!
567 !! FUNCTION
568 !!
569 !! Initialize population energy and gradient arrays to zero before any
570 !! possible MPI_Allreduce calls are made
571 !!
572 !! SYNOPSIS
573 subroutine dlf_init_engarrays
574 !! SOURCE
576 ! This (re)initialization is necessary for the parallel implementation, when
577 ! MPI_Allreduce calls are made.
578 if (stochastic) then
579 pop_xgradient(:,:,:) = 0.0D0
580 pop_energies(:) = 0.0D0
581 else if (genetic) then
582 ! the first member of the population is excluded from mutation and is
583 ! thus unchanged in this generation - no need to recalculate
584 if (glob%mytask == 0) then
585 pop_xgradient(2:glob%po_pop_size,:,:) = 0.0D0
586 pop_energies(2:glob%po_pop_size) = 0.0D0
587 else
588 pop_xgradient(:,:,:) = 0.0D0
589 pop_energies(:) = 0.0D0
590 end if
591 end if
593 end subroutine dlf_init_engarrays
594 !!****
597 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
598 !!****f* parallel_opt/dlf_get_engarrays
599 !!
600 !! FUNCTION
601 !!
602 !! Calculate energy and gradient for each individual in the population
603 !! (apart from the first individual in a GA run...)
604 !! Parallelised, with each taskfarm essentially calculating e and g for
605 !! the i-th, (i+ntasks)-th, (i+2*ntasks)-th, etc. individuals.
606 !!
607 !! INPUTS
608 !!
609 !! pop_icoords(:,:)
610 !! lower_index
611 !! glob%po_pop_size
612 !! glob%ntasks
613 !! glob%mytask
614 !! iimage
615 !! glob%nvar
616 !!
617 !! OUTPUTS
618 !!
619 !! pop_energies(:)
620 !! pop_xgradient(:,:,:)
621 !!
622 !! SYNOPSIS
623 subroutine dlf_get_engarrays
624 !! SOURCE
626 integer :: l, m, status
628 do m = lower_index, glob%po_pop_size
629 l = mod(m,glob%ntasks)
630 stat%sene = stat%sene + 1
631 if (l == glob%mytask) then
633 glob%icoords(:) = pop_icoords(m,:)
635 ! convert glob%icoords (used as a temporary storage array) back
636 ! to Cartesians in glob%xcoords, before call to dlf_get_gradient
638 call clock_start("COORDS")
639 call dlf_coords_itox(iimage)
640 call clock_stop("COORDS")
642 call clock_start("EANDG")
643 call dlf_get_gradient(glob%nvar,glob%xcoords,pop_energies(m), &
644 pop_xgradient(m,:,:),iimage,&
645 #ifdef GAMESS
646 core,&
647 #endif
648 status)
649 stat%pene = stat%pene + 1
650 call clock_stop("EANDG")
652 ! check of NaN in the energy (comparison of NaN with any number
653 ! returns .false. , pgf90 does not understand isnan() )
654 if ( abs(pop_energies(m)) > huge(1.D0) ) then
655 status = 1
656 else
657 if (.not. abs(pop_energies(m)) < huge(1.D0) ) status = 1
658 end if
660 if (status/=0) then
661 pop_energies(m) = huge(1.0D0)
662 pop_xgradient(m,:,:) = huge(1.0D0)
663 end if
665 endif
666 enddo
668 end subroutine dlf_get_engarrays
669 !!****
672 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
673 !!****f* parallel_opt/dlf_stoch_makepop
674 !!
675 !! FUNCTION
676 !!
677 !! Create the population for a stochastic search.
678 !! Only executed on the rank-zero processor, as it depends on random
679 !! numbers.
680 !! Works wholly in "i" coordinates.
681 !!
682 !! INPUTS
683 !!
684 !! glob%po_pop_size
685 !! glob%nivar
686 !! glob%iam
687 !! glob%po_distribution
688 !! glob%po_radius(:)
689 !! glob%po_scalefac
690 !! igradient_best(:)
691 !!
692 !! OUTPUTS
693 !!
694 !! pop_icoords(:,:)
695 !!
696 !! SYNOPSIS
697 subroutine dlf_stoch_makepop
698 !! SOURCE
700 integer :: l,m
702 if (glob%iam /= 0) return ! as only iam == 0 has seeded the random number
703 ! generator.
705 ! we are working wholly in "i" coordinates in this subroutine
707 pop_icoords(:,:) = 0.0D0
708 do l = 1, glob%po_pop_size
709 do m = 1, glob%nivar
710 call random_number(random)
711 select case (glob%po_distribution)
712 case (3) ! force bias
713 random = random * glob%po_radius(m)
714 random = random * glob%po_scalefac * abs(igradient_best(m))
715 if (igradient_best(m) > 0.0D0) random = -1.0D0*random
716 case (2) ! force_direction_bias
717 random = random * glob%po_radius(m)
718 if (igradient_best(m) > 0.0D0) random = -1.0D0*random
719 case (1)
720 ! uniform distribution over the full hypersphere in config.space
721 random = (random * 2.0D0 * glob%po_radius(m)) - glob%po_radius(m)
722 case default
723 write(stderr,'(1x,a,i2,a)') "Parallel optimisation &
724 &distribution type ",glob%po_distribution," not implemented"
725 call dlf_fail("Bad option for parallel optimisation distribution &
726 &type")
727 end select
728 pop_icoords(l,m) = icoords_best(m) + random
729 end do
730 enddo
732 end subroutine dlf_stoch_makepop
733 !!****
736 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
737 !!****f* parallel_opt/dlf_genetic_makepop
738 !!
739 !! FUNCTION
740 !!
741 !! Create the population during the cycles of the genetic algorithm.
742 !! Only executed on the rank-zero processor, as it depends on random
743 !! numbers.
744 !! Works wholly in "i" coordinates.
745 !!
746 !! INPUTS
747 !!
748 !! glob%iam
749 !! glob%nivar
750 !! glob%po_pop_size
751 !! noffspring
752 !! pop_energies(:)
753 !! pop_icoords(:,:)
754 !!
755 !! OUTPUTS
756 !!
757 !! tno_diversity
758 !! pop_icoords(:,:)
759 !!
760 !! SYNOPSIS
761 subroutine dlf_genetic_makepop
762 !! SOURCE
764 integer :: l, m, kk, cross_over
765 integer,dimension(2) :: ip
766 real(rk) :: sum_w, dummy, dummy2, delta, harvest, blending
767 real(rk), allocatable :: w_cost(:)
769 if (glob%iam /= 0) return ! as only iam == 0 has seeded the random number
770 ! generator.
772 ! set cost(i.e. energy)-weighted probabilities of parenthood for the breeding
773 ! population, i.e. the glob%po_pop_size - noffspring lowest-energy individuals
775 call allocate(w_cost,glob%po_pop_size - noffspring)
776 sum_w = 0.0D0
778 do l = 1, glob%po_pop_size - noffspring
779 w_cost(l) = pop_energies(l) - pop_energies(glob%po_pop_size - noffspring + 1)
780 end do
782 sum_w = sum(w_cost)
784 ! test whether the (sorted) energies of the breeding section of the
785 ! population are all identical to at least a chosen precision.
787 dummy = abs( (pop_energies(1) - pop_energies(glob%po_pop_size-noffspring)) &
788 & / pop_energies(1) )
790 ! Also test for the case where, of the breeding population plus the
791 ! next-highest-energy configuration, the lowest is unique and the rest are
792 ! identical. This situation would mean that the lowest-energy individual
793 ! is the first parent, and no second parent can be selected in the scheme
794 ! in the offspring generation loop below.
796 dummy2 = abs( (pop_energies(2) - pop_energies(glob%po_pop_size-noffspring+1)) &
797 & / pop_energies(2) )
799 ! if the population is lacking in genetic diversity then return
801 if (sum_w == 0.0D0 .or. dummy < small .or. dummy2 < small) then
802 tno_diversity = .true.
803 call deallocate(w_cost)
804 return
805 end if
807 tno_diversity = .false.
809 ! Accumulate the w_cost's
811 w_cost(:) = w_cost(:)/sum_w
813 do l = 2, glob%po_pop_size - noffspring
814 w_cost(l) = w_cost(l) + w_cost(l-1)
815 end do
817 ! in case there are rounding errors, set the final w_cost to 1.0
819 w_cost(glob%po_pop_size - noffspring) = 1.0D0
821 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
822 ! offspring generation loop
823 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
825 do l = 1, noffspring - 1, 2
827 ! select parents: ip(1:2). Makes sure that ip(1) /= ip(2)
829 kk = 1
830 do
831 call random_number(harvest)
832 if (harvest < w_cost(1)) then
833 ip(kk) = 1
834 kk = kk + 1
835 else
836 do m = 2, glob%po_pop_size - noffspring
837 if (harvest <= w_cost(m)) then
838 ip(kk) = m
839 kk = kk + 1
840 exit
841 end if
842 end do
843 end if
844 if (kk == 3) then
845 if (ip(1) == ip(2)) then
846 kk = 2 ! choose a different "father"
847 else
848 exit ! we have two different parents
849 end if
850 end if
851 end do
853 ! choose one of the glob%nivar variables as the cross-over point
855 call random_number(harvest)
856 cross_over = int( harvest * glob%nivar ) + 1
857 delta = pop_icoords(ip(1),cross_over) - pop_icoords(ip(2),cross_over)
859 ! Check we have a viable cross over point given the parents
861 if ( (abs( delta/pop_icoords(ip(1),cross_over) ) < 1.0D3*small) .and. &
862 &(cross_over == 1 .or. cross_over == glob%nivar) ) then
863 ! no point crossing as the offspring will be identical to the parents
865 ! Shift the cross-over point one element to the right until we find a
866 ! point at which the parents differ.
867 ! At present, if no such point exists then we just continue with the
868 ! value of cross_over from after the do loop, regardless. IMPROVE !!!
869 do m = 1, glob%nivar - 1
870 cross_over = cross_over + 1
871 if (cross_over > glob%nivar) cross_over = 1
872 delta = pop_icoords(ip(1),cross_over)-pop_icoords(ip(2),cross_over)
873 if (abs( delta/pop_icoords(ip(1),cross_over) ) > 1.0D3*small) exit
874 end do
875 end if
877 ! Choose blending factor
879 call random_number(blending)
881 ! Make offspring
883 m = glob%po_pop_size - noffspring ! to save characters
885 pop_icoords(m + l , cross_over) = pop_icoords(ip(1), cross_over) &
886 & - blending * delta
887 pop_icoords(m + l + 1 , cross_over) = pop_icoords(ip(2), cross_over) &
888 & + blending * delta
889 if (cross_over == glob%nivar) then
890 pop_icoords(m + l , 1 : glob%nivar - 1) = &
891 & pop_icoords(ip(1), 1 : glob%nivar - 1)
892 pop_icoords(m + l + 1 , 1 : glob%nivar - 1) = &
893 & pop_icoords(ip(2), 1 : glob%nivar - 1)
894 else
895 if (cross_over == 1) then
896 pop_icoords(m + l , 2 : glob%nivar) = &
897 & pop_icoords(ip(2), 2 : glob%nivar)
898 pop_icoords(m + l + 1 , 2 : glob%nivar) = &
899 & pop_icoords(ip(1), 2 : glob%nivar)
900 else
901 pop_icoords(m + l , 1 : cross_over - 1) = &
902 & pop_icoords(ip(1), 1 : cross_over - 1)
903 pop_icoords(m + l + 1 , 1 : cross_over - 1) = &
904 & pop_icoords(ip(2), 1 : cross_over - 1)
906 pop_icoords(m + l , cross_over + 1 : glob%nivar) = &
907 & pop_icoords(ip(2), cross_over + 1 : glob%nivar)
908 pop_icoords(m + l + 1 , cross_over + 1 : glob%nivar) = &
909 & pop_icoords(ip(1), cross_over + 1 : glob%nivar)
910 end if
911 end if
912 end do
914 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
915 ! end of offspring generation loop
916 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
918 call deallocate(w_cost)
920 end subroutine dlf_genetic_makepop
921 !!****
924 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
925 !!****f* parallel_opt/dlf_genetic_initialpop
926 !!
927 !! FUNCTION
928 !!
929 !! Create the "initial" population for the genetic algorithm, including
930 !! the population after a reset has been requested, from a starting
931 !! structure in icoords form.
932 !! Only executed on the rank-zero processor, as it depends on random
933 !! numbers.
934 !!
935 !! INPUTS
936 !!
937 !! princeps(:)
938 !! n
939 !! glob%nivar
940 !! glob%po_radius(:)
941 !! glob%iam
942 !!
943 !! OUTPUTS
944 !!
945 !! dum_pop(:,:)
946 !!
947 !! SYNOPSIS
948 subroutine dlf_genetic_initialpop(dum_pop,n,princeps)
949 !! SOURCE
951 ! princeps is the "first citizen" in icoords form, from which the rest of the
952 ! population will be seeded.
953 ! n is the population size (can be either the initial or working population)
954 ! On return, dum_pop contains the newly created population in icoords form
956 integer :: n, l, m
957 real(rk) :: r
958 real(rk) :: princeps(:), dum_pop(:,:)
960 if (glob%iam /= 0) return ! as only iam == 0 has seeded the random number
961 ! generator.
963 dum_pop(1,:) = princeps(:)
965 do l = 2, n ! n is the population size
966 do m = 1, glob%nivar
967 call random_number(r)
968 r = (r * 2.0D0 * glob%po_radius(m)) - glob%po_radius(m)
969 dum_pop(l,m) = dum_pop(1,m) + r
970 end do
971 end do
973 end subroutine dlf_genetic_initialpop
974 !!****
977 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
978 !!****f* parallel_opt/dlf_genetic_mutate
979 !!
980 !! FUNCTION
981 !!
982 !! Mutate the population in the genetic algorithm by perturbing i-coordinates.
983 !! Only executed on the rank-zero processor, as it depends on random
984 !! numbers.
985 !!
986 !! INPUTS
987 !!
988 !! glob%iam
989 !! nmutations
990 !! glob%po_pop_size
991 !! glob%nivar
992 !! glob%po_radius(:)
993 !! pop_icoords(:,:)
994 !!
995 !! OUTPUTS
996 !!
997 !! pop_icoords(:,:)
998 !!
999 !! SYNOPSIS
1000 subroutine dlf_genetic_mutate
1001 !! SOURCE
1003 integer :: l, indiv, element
1004 real(rk) :: harvest, shift, factor
1006 if (glob%iam /= 0) return ! as only iam == 0 has seeded the random number
1007 ! generator.
1009 do l = 1, nmutations
1010 call random_number(harvest)
1011 indiv = int( harvest*(glob%po_pop_size - 1) ) + 2
1012 ! the "+ 2" means indiv is never one (i.e. the lowest-energy individual)
1013 if (indiv > glob%po_pop_size) indiv = glob%po_pop_size
1015 call random_number(harvest)
1016 element = int( harvest*glob%nivar ) + 1
1017 if (element > glob%nivar) element = glob%nivar
1019 call random_number(harvest)
1020 factor = (glob%po_radius(element)*glob%po_radius(element))/2.0D0
1021 shift = -factor * log( sqrt( 3.1415D0*factor )*harvest )
1022 !!! why this form?
1023 if (abs(shift) > glob%po_radius(element)) shift = glob%po_radius(element)
1025 call random_number(harvest)
1026 if (harvest < 0.5D0) shift = -shift
1027 pop_icoords(indiv, element) = pop_icoords(indiv, element) + shift
1028 end do
1030 end subroutine dlf_genetic_mutate
1031 !!****
1034 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1035 !!****f* parallel_opt/dlf_stoch_endcycle
1036 !!
1037 !! FUNCTION
1038 !!
1039 !! Book-keeping at the end of a stochastic search cycle, including
1040 !! tests for failure and convergence and checkpointing.
1041 !!
1042 !! INPUTS
1043 !!
1044 !! pop_energies(:)
1045 !! pop_icoords(:,:)
1046 !! pop_xgradient(:,:,:)
1047 !! glob%po_tolerance_g
1048 !! glob%po_tolerance_r(:)
1049 !! glob%iam
1050 !! glob%po_radius(:)
1051 !! glob%po_contraction
1052 !! glob%dump
1053 !! printl
1054 !! printf
1055 !!
1056 !! OUTPUTS
1057 !!
1058 !! energy_best
1059 !! icoords_best(:)
1060 !! xcoords_best(:,:)
1061 !! igradient_best(:)
1062 !! tconv
1063 !! texit
1064 !!
1065 !! files, if converged -- best.xyz
1066 !! best_active.xyz
1067 !!
1068 !! SYNOPSIS
1069 subroutine dlf_stoch_endcycle
1070 !! SOURCE
1072 integer :: l, k
1074 if ( minval(pop_energies(:)) < energy_best ) then ! store new "best" arrays
1075 energy_best = minval(pop_energies(:))
1076 position(1:1) = minloc(pop_energies(:))
1077 icoords_best(:) = pop_icoords(position(1),:)
1079 ! convert icoords_best(:) to Cartesians and store in xcoords_best(3,:)
1080 glob%icoords(:) = icoords_best(:)
1081 call clock_start("COORDS")
1082 call dlf_coords_itox(iimage)
1083 call clock_stop("COORDS")
1084 xcoords_best(:,:) = glob%xcoords(:,:)
1086 ! convert pop_xgradient(position(1), :, :) from Cartesians and store in
1087 ! igradient_best(:)
1088 glob%xgradient(:,:) = pop_xgradient(position(1), :, :)
1089 call clock_start("COORDS")
1090 call dlf_coords_xtoi(dummy_logical,dummy_logical,iimage,0)
1091 call clock_stop("COORDS")
1092 igradient_best(:) = glob%igradient(:)
1093 endif
1095 write(stdout,'(1x,a,i10)')"CYCLE ",i
1096 ! i is the counter in the cycles do loop of the calling subroutine
1098 ! Make sure the value of glob%po_tolerance_g is appropriate for the coord system
1099 if ( maxval(abs(igradient_best(:))) < glob%po_tolerance_g ) then
1100 if (printl > 0) then
1101 write(stdout,'(1x,a)')"****Stochastic search converged****"
1102 write(stdout,'(1x,a,es16.9)')"Lowest energy = ",energy_best
1103 if (printl >= 4) then
1104 write(stdout,'(1x,a)')"for Cartesian coordinates "
1105 write(stdout,'(3(1x,es16.9))') xcoords_best(:,:)
1106 end if
1107 write(stdout,'(1x,a,es16.9)')"Max component of absolute force = ",&
1108 &maxval(abs(igradient_best(:)))
1109 write(stdout,'(1x,a,i10)')"Number of energy evaluations required = ",&
1110 &1 + glob%po_pop_size*i
1111 end if
1113 if (printf>=2 .and. glob%iam == 0) then
1114 call clock_start("XYZ")
1115 open(unit=400,file="best.xyz")
1116 open(unit=401,file="best_active.xyz")
1117 call write_xyz(400,glob%nat,glob%znuc,xcoords_best)
1118 call write_xyz_active(401,glob%nat,glob%znuc,glob%spec,xcoords_best)
1119 close(400)
1120 close(401)
1121 call dlf_put_coords(glob%nvar,1,energy_best,xcoords_best,glob%iam)
1122 call clock_stop("XYZ")
1123 end if
1124 tconv = .true.
1125 texit = .true.
1126 return
1127 end if
1129 glob%po_radius(:) = glob%po_contraction * glob%po_radius(:)
1131 ! Stop if any component of the vector of radii falls below its tolerance.
1132 do l = 1, glob%nivar
1133 if (glob%po_radius(l) < glob%po_tolerance_r(l)) then
1134 if (printl > 0) then
1135 write(stdout,'(1x,a)')"Search restricted to less than minimum radius"
1136 write(stdout,'(1x,a)')"Sample radius, tolerance:"
1137 do k = 1, glob%nivar
1138 write(stdout,'(2(1x,es16.9))') glob%po_radius(k),glob%po_tolerance_r(k)
1139 end do
1140 write(stdout,'(1x,a,es16.9)')"Lowest energy = ",energy_best
1141 if (printl >= 4) then
1142 write(stdout,'(1x,a)')"for Cartesian coordinates "
1143 write(stdout,'(3(1x,es16.9))') xcoords_best(:,:)
1144 end if
1145 write(stdout,'(1x,a,es16.9)')"Max component of absolute force = ", &
1146 &maxval(abs(igradient_best(:)))
1147 write(stdout,'(1x,a,i10)')"Number of energy evaluations required = ",&
1148 &1 + glob%po_pop_size*i
1149 write(stdout,'(1x,a)')"Stopping"
1150 end if
1151 texit = .true.
1152 return
1153 end if
1154 end do
1156 if (glob%dump > 0 .and. mod(i,glob%dump)==0) then
1157 call clock_start("CHECKPOINT")
1158 if (printl>=6) write(stdout,"('Writing restart information')")
1159 status = 1
1160 call dlf_checkpoint_write(status)
1161 call dlf_checkpoint_po_write
1162 call clock_stop("CHECKPOINT")
1163 end if
1165 write(stdout,'(1x,a,es16.9)')"Lowest energy = ",energy_best
1166 if (printl >= 4) then
1167 write(stdout,'(1x,a)')"for Cartesian coordinates "
1168 write(stdout,'(3(1x,es16.9))') xcoords_best(:,:)
1169 end if
1170 write(stdout,'(1x,a,es16.9)')"Max component of absolute force = ", &
1171 &maxval(abs(igradient_best(:)))
1172 write(stdout,'(1x,a,i10)')"Energy calls = ",1 + glob%po_pop_size*i
1174 call flush(stdout)
1176 end subroutine dlf_stoch_endcycle
1177 !!****
1180 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1181 !!****f* parallel_opt/dlf_genetic_endcycle
1182 !!
1183 !! FUNCTION
1184 !!
1185 !! Book-keeping at the end of a genetic algorithm cycle, including
1186 !! checkpointing and a test for convergence.
1187 !!
1188 !! INPUTS
1189 !!
1190 !! pop_energies(:)
1191 !! pop_icoords(:,:)
1192 !! pop_xgradient(:,:,:)
1193 !! glob%po_tolerance_g
1194 !! glob%iam
1195 !! glob%po_nsave
1196 !! glob%dump
1197 !! printl
1198 !! printf
1199 !!
1200 !! OUTPUTS
1201 !!
1202 !! energy_best
1203 !! icoords_best(:)
1204 !! xcoords_best(:,:)
1205 !! igradient_best(:)
1206 !! tconv
1207 !!
1208 !! if nsave/=0 -- energies_save(:)
1209 !! nevals_save(:)
1210 !! xcoords_save(:,:,:)
1211 !!
1212 !! local variables -- mean_e
1213 !! sigma_e
1214 !!
1215 !! files, if converged -- best.xyz
1216 !! best_active.xyz
1217 !!
1218 !! SYNOPSIS
1219 subroutine dlf_genetic_endcycle
1220 !! SOURCE
1222 integer :: m, l
1223 real(rk) :: mean_e, mean_e2, sigma_e
1225 energy_best = pop_energies(1)
1226 icoords_best(:) = pop_icoords(1,:)
1228 ! convert icoords_best(:) to Cartesians and store in xcoords_best(3,:)
1229 glob%icoords(:) = icoords_best(:)
1230 call clock_start("COORDS")
1231 call dlf_coords_itox(iimage)
1232 call clock_stop("COORDS")
1233 xcoords_best(:,:) = glob%xcoords(:,:)
1235 ! convert pop_xgradient(1, :, :) from Cartesians and store in
1236 ! igradient_best(:) !!! assuming we know the gradient already....
1237 glob%xgradient(:,:) = pop_xgradient(1, :, :)
1238 call clock_start("COORDS")
1239 call dlf_coords_xtoi(dummy_logical,dummy_logical,iimage,0)
1240 call clock_stop("COORDS")
1241 igradient_best(:) = glob%igradient(:)
1243 mean_e = sum(pop_energies(:)) / dble(glob%po_pop_size)
1244 mean_e2 = sum(pop_energies(:)**2) / dble(glob%po_pop_size)
1245 sigma_e = sqrt(mean_e2 - mean_e**2)
1247 ! Write details for this generation
1248 ! i is the counter in the cycles do loop of the calling subroutine
1249 write(stdout,'(1x,a,i10)')"CYCLE ",i
1250 write(stdout,'(1x,a,es16.9)')"Lowest energy = ",pop_energies(1)
1251 if (printl >= 4) then
1252 write(stdout,'(1x,a)')"for Cartesian coordinates "
1253 write(stdout,'(3(1x,es16.9))') xcoords_best(:,:)
1254 end if
1255 write(stdout,'(1x,a,es16.9)')"Max component of absolute force = ",&
1256 &maxval(abs(igradient_best(:)))
1257 write(stdout,'(1x,a,es16.9)')"Mean energy for population = ",mean_e
1258 write(stdout,'(1x,a,es16.9)')"Standard deviation of the energy = ",sigma_e
1259 write(stdout,'(1x,a,i10)')"Energy calls = ", &
1260 &glob%po_init_pop_size + (glob%po_pop_size - 1)*(i + nresets)
1262 ! Make sure the value of glob%po_tolerance_g is appropriate for the coord system
1263 if ( (maxval(abs(igradient_best(:))) < glob%po_tolerance_g) .AND. &
1264 (energy_best < HUGE(1.0D0)) ) then
1265 ! Added a catch for energy_best == huge above to avoid false convergence in the (unlikely!) case the coord transformation
1266 ! for the gradient vastly changes the values from when the Cartesian gradient vector was set to huge in dlf_get_engarrays
1267 if (printl > 0) then
1268 write(stdout,'(1x,a)')"****Genetic algorithm converged to a minimum****"
1269 write(stdout,'(1x,a,es16.9)')"Lowest energy = ",energy_best
1270 if (printl >= 4) then
1271 write(stdout,'(1x,a)')"for Cartesian coordinates "
1272 write(stdout,'(3(1x,es16.9))') xcoords_best(:,:)
1273 end if
1274 write(stdout,'(1x,a,es16.9)')"Max component of absolute force = ", &
1275 &maxval(abs(igradient_best(:)))
1276 write(stdout,'(1x,a,i10)')"Number of energy evaluations required = ",&
1277 &glob%po_init_pop_size + (glob%po_pop_size - 1)*(i + nresets)
1278 write(stdout,'(1x,a)')"GA continuing... (minimum may be only local)"
1279 end if
1281 if (printf>=2 .and. glob%iam == 0) then
1282 call clock_start("XYZ")
1283 open(unit=400,file="best.xyz",position="APPEND")
1284 open(unit=401,file="best_active.xyz",position="APPEND")
1285 call write_xyz(400,glob%nat,glob%znuc,xcoords_best)
1286 call write_xyz_active(401,glob%nat,glob%znuc,glob%spec,xcoords_best)
1287 close(400)
1288 close(401)
1289 call clock_stop("XYZ")
1290 end if
1292 if (glob%po_nsave > 0) then
1293 if (energy_best < energies_save(glob%po_nsave)) then
1294 do m = 1, glob%po_nsave
1295 !!!crudely check for multiple incidences of a particular minimum
1296 if ( abs( ( energy_best - energies_save(m) )/energy_best ) &
1297 & < small ) exit
1298 if (energy_best < energies_save(m)) then
1299 do l = glob%po_nsave-1, m, -1
1300 energies_save(l+1) = energies_save(l)
1301 xcoords_save(l+1,:,:) = xcoords_save(l,:,:)
1302 nevals_save(l+1) = nevals_save(l)
1303 end do
1304 energies_save(m) = energy_best
1305 xcoords_save(m,:,:) = xcoords_best(:,:)
1306 nevals_save(m) = glob%po_init_pop_size + &
1307 &(glob%po_pop_size - 1)*(i + nresets)
1308 exit
1309 end if
1310 end do
1311 end if
1312 end if
1314 tconv = .true.
1315 end if
1317 if (glob%dump > 0 .and. mod(i,glob%dump)==0) then
1318 call clock_start("CHECKPOINT")
1319 if (printl>=6) write(stdout,"('Writing restart information')")
1320 status = 1
1321 call dlf_checkpoint_write(status)
1322 call dlf_checkpoint_po_write
1323 call clock_stop("CHECKPOINT")
1324 end if
1326 call flush(stdout)
1328 end subroutine dlf_genetic_endcycle
1329 !!****
1332 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1333 !!****f* parallel_opt/dlf_po_destroy
1334 !!
1335 !! FUNCTION
1336 !!
1337 !! Deallocate population arrays
1338 !!
1339 !! SYNOPSIS
1340 subroutine dlf_po_destroy
1341 !! SOURCE
1343 call deallocate(pop_icoords)
1344 call deallocate(pop_xgradient)
1345 call deallocate(pop_energies)
1347 call deallocate(icoords_best)
1348 call deallocate(igradient_best)
1349 call deallocate(xcoords_best)
1351 if (genetic) then
1353 call deallocate(init_pop_icoords)
1354 call deallocate(init_pop_energies)
1356 end if
1358 call deallocate(glob%po_radius)
1359 call deallocate(glob%po_tolerance_r)
1361 end subroutine dlf_po_destroy
1362 !!****
1365 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1366 !!****f* parallel_opt/dlf_checkpoint_po_read
1367 !!
1368 !! FUNCTION
1369 !!
1370 !! Reading of the parallel optimisation checkpoint file
1371 !!
1372 !! OUTPUTS
1373 !!
1374 !! trestarted
1375 !!
1376 !! SS -- energy_best
1377 !! icoords_best(:)
1378 !! xcoords_best(:,:)
1379 !! igradient_best(:)
1380 !!
1381 !! GA -- pop_icoords(:,:)
1382 !! pop_energies(:)
1383 !! pop_xgradient(:,:,:)
1384 !!
1385 !! SYNOPSIS
1386 subroutine dlf_checkpoint_po_read
1387 !! SOURCE
1389 logical :: tchk, tallocated
1390 real(rk), allocatable :: dummy(:,:,:)
1392 trestarted=.false.
1394 ! check if checkpoint file exists
1395 INQUIRE(FILE="dlf_parallel_opt.chk",EXIST=tchk)
1396 if (.not.tchk) then
1397 write(stdout,10) "File dlf_parallel_opt.chk not found"
1398 return
1399 end if
1401 ! open the checkpoint file
1402 if (tchkform) then
1403 open(unit=100,file="dlf_parallel_opt.chk",form="formatted")
1404 else
1405 open(unit=100,file="dlf_parallel_opt.chk",form="unformatted")
1406 end if
1408 if (stochastic) then
1410 call read_separator(100,"Best arrays",tchk)
1411 if (.not. tchk) return
1412 if (tchkform) then
1413 read(100,*,end=201,err=200) energy_best, icoords_best(:), &
1414 &xcoords_best(:,:), igradient_best(:)
1415 else
1416 read(100,end=201,err=200) energy_best, icoords_best(:), &
1417 &xcoords_best(:,:), igradient_best(:)
1418 end if
1420 else if (genetic) then
1422 call read_separator(100,"Population",tchk)
1423 if (.not. tchk) return
1424 if (tchkform) then
1425 read(100,*,end=201,err=200) pop_icoords(:,:), pop_energies(:)
1426 read(100,*,end=201,err=200) tallocated
1427 if (tallocated .and. allocated(pop_xgradient)) then
1428 read(100,*,end=201,err=200) pop_xgradient(:,:,:)
1429 else if (tallocated) then
1430 call allocate(dummy, glob%po_pop_size, 3, glob%nat)
1431 read(100,*,end=201,err=200) dummy(:,:,:)
1432 call deallocate(dummy)
1433 end if
1434 else
1435 read(100,end=201,err=200) pop_icoords(:,:), pop_energies(:)
1436 read(100,end=201,err=200) tallocated
1437 if (tallocated .and. allocated(pop_xgradient)) then
1438 read(100,end=201,err=200) pop_xgradient(:,:,:)
1439 else if (tallocated) then
1440 call allocate(dummy, glob%po_pop_size, 3, glob%nat)
1441 read(100,end=201,err=200) dummy(:,:,:)
1442 call deallocate(dummy)
1443 end if
1444 end if
1446 end if
1448 call read_separator(100,"END",tchk)
1449 if (.not.tchk) return
1451 ! successful reading
1452 close(100)
1453 trestarted=.true.
1454 return
1456 ! return on error
1457 200 continue
1458 write(stdout,10) "Error reading parallel optimisation checkpoint file"
1459 return
1460 201 continue
1461 write(stdout,10) "Error (EOF) reading parallel optimisation checkpoint file"
1462 return
1464 10 format("Checkpoint reading WARNING: ",a)
1466 end subroutine dlf_checkpoint_po_read
1467 !!****
1470 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1471 !!****f* parallel_opt/dlf_checkpoint_po_write
1472 !!
1473 !! FUNCTION
1474 !!
1475 !! Writing of the parallel optimisation checkpoint file
1476 !!
1477 !! INPUTS
1478 !!
1479 !! trestarted
1480 !!
1481 !! SS -- energy_best
1482 !! icoords_best(:)
1483 !! xcoords_best(:,:)
1484 !! igradient_best(:)
1485 !!
1486 !! GA -- pop_icoords(:,:)
1487 !! pop_energies(:)
1488 !! pop_xgradient(:,:,:)
1489 !!
1490 !! SYNOPSIS
1491 subroutine dlf_checkpoint_po_write
1492 !! SOURCE
1494 ! Only want one processor to do the writing; that processor must know
1495 ! all the necessary information.
1496 if (glob%iam /= 0) return
1498 ! Open the checkpoint file
1499 if (tchkform) then
1500 open(unit=100,file="dlf_parallel_opt.chk",form="formatted")
1501 else
1502 open(unit=100,file="dlf_parallel_opt.chk",form="unformatted")
1503 end if
1505 if (stochastic) then
1507 call write_separator(100,"Best arrays")
1508 if (tchkform) then
1509 write(100,*) energy_best, icoords_best(:), &
1510 &xcoords_best(:,:), igradient_best(:)
1511 else
1512 write(100) energy_best, icoords_best(:), &
1513 &xcoords_best(:,:), igradient_best(:)
1514 end if
1516 else if (genetic) then
1518 call write_separator(100,"Population")
1519 if (tchkform) then
1520 write(100,*) pop_icoords(:,:), pop_energies(:)
1521 write(100,*) allocated(pop_xgradient)
1522 if (allocated(pop_xgradient)) write(100,*) pop_xgradient(:,:,:)
1523 else
1524 write(100) pop_icoords(:,:), pop_energies(:)
1525 write(100) allocated(pop_xgradient)
1526 if (allocated(pop_xgradient)) write(100) pop_xgradient(:,:,:)
1527 end if
1529 end if
1531 call write_separator(100,"END")
1532 close(100)
1534 end subroutine dlf_checkpoint_po_write
1535 !!****
1538 end subroutine dlf_parallel_opt
1539 !!****
