Line data Source code
1 : !!****m* ABINIT/predict_ga
2 : !! NAME
3 : !! predict_ga
4 : !!
5 : !! FUNCTION
6 : !! Given a given set of images, which represent a population, it predicts a new set of images.
7 : !! The implementation is based on a Genetic Algorithm idea, where the best fit candidates are passed
8 : !! to the next generation. Those are chosen from ga_opt_percent% best fit and (1-ga_opt_percent)% from Genetic rules
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2009-2026 ABINIT group (XG, AHR)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : !!
17 : !! INPUTS
18 : !! itimimage=time index for image propagation (itimimage+1 is to be predicted here)
19 : !! itimimage_eff=time index in the history
20 : !! list_dynimage(nimage)=list of dynamical images.
21 : !! This is quite useful when ground states of the A and B states is known
22 : !! natom=dimension of vel_timimage and xred_timimage
23 : !! ndynimage=number of dynamical images
24 : !! nimage= population size
25 : !! ntimimage_stored=number of time steps stored in the history
26 : !!
27 : !! OUTPUT
28 : !!
29 : !! SIDE EFFECTS
30 : !! results_img(ntimimage_stored,nimage)=datastructure that holds the history of previous computations.
31 : !! results_img(:,:)%acell(3)
32 : !! at input, history of the values of acell for all images
33 : !! at output, the predicted values of acell for all images
34 : !! results_img(:,:)%results_gs
35 : !! at input, history of the values of energies and forces for all images
36 : !! results_img(:,:)%rprim(3,3)
37 : !! at input, history of the values of rprim for all images
38 : !! at output, the predicted values of rprim for all images
39 : !! results_img(:,:)%vel(3,natom)
40 : !! at input, history of the values of vel for all images
41 : !! at output, the predicted values of vel for all images
42 : !! results_img(:,:)%vel_cell(3,3)
43 : !! at input, history of the values of vel_cell for all images
44 : !! at output, the predicted values of vel_cell for all images
45 : !! results_img(:,:)%xred(3,natom)
46 : !! at input, history of the values of xred for all images
47 : !! at output, the predicted values of xred for all images
48 : !!
49 : !! SOURCE
50 :
51 : #if defined HAVE_CONFIG_H
52 : #include "config.h"
53 : #endif
54 :
55 : #include "abi_common.h"
56 :
57 : MODULE m_use_ga
58 :
59 : use defs_basis
60 : use m_abicore
61 : use m_ga
62 : use m_sort
63 :
64 : use m_symfind, only : symfind, symanal, symlatt
65 : use m_geometry, only : mkradim, mkrdim, metric, dist2
66 : use m_results_img, only : results_img_type,gather_array_img
67 : use m_numeric_tools, only : uniformrandom
68 :
69 : implicit none
70 :
71 : private
72 :
73 : public :: predict_ga
74 :
75 : CONTAINS
76 :
77 4 : subroutine predict_ga(itimimage_eff,idum,ga_param,natom,nimage,&
78 4 : & ntimimage_stored,results_img)
79 :
80 : !Arguments ------------------------------------
81 : !scalars
82 : integer,intent(in) :: itimimage_eff,natom,nimage,ntimimage_stored
83 : integer,intent(inout) :: idum
84 : !arrays
85 : type(results_img_type) :: results_img(nimage,ntimimage_stored)
86 : type(ga_type),intent(inout) :: ga_param
87 :
88 : !Local variables------------------------------
89 : !scalars
90 : integer :: ii,jj,kk,itmp,iimage,indiv,oper,iparent1,iparent2,ndimen
91 : integer :: next_itimimage,nsurvivor,nspinor
92 : !character(len=500) :: message
93 : !arrays
94 4 : integer,allocatable :: ieperm(:),ihperm(:),ibgoptperm(:,:),ibgperm(:,:)
95 4 : integer,allocatable :: zperm1(:),zperm2(:)
96 : real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3)
97 : real(dp) :: rprimdparent1(3,3),rprimdparent2(3,3)
98 : real(dp) :: rprimson1(3,3),rprimson2(3,3)
99 : real(dp) :: rprimdson1(3,3),rprimdson2(3,3)
100 : real(dp) :: acellson1(3),acellson2(3)
101 8 : real(dp) :: son1(3,natom),son2(3,natom)
102 8 : real(dp) :: parent1(3,natom),parent2(3,natom)
103 : ! store the energy and the enthalpy of all elements of the population
104 4 : real(dp),allocatable :: etotal_img(:),enthalpy_img(:),bg_img(:,:),bg_opt_img(:,:)
105 4 : real(dp),allocatable :: acell(:,:),acell_old(:,:),rprim(:,:,:),rprim_old(:,:,:),rprimd(:,:,:)
106 4 : real(dp),allocatable :: fitness(:),zcoor1(:),zcoor2(:)
107 4 : real(dp),allocatable :: coor(:,:),coor_old(:,:)
108 4 : real(dp),allocatable :: vson1(:),vson2(:)
109 :
110 : !real quantities
111 : real(dp) :: denom,sumH,Hmin,Hmax,ucvol,rtmp(3,3)
112 :
113 : ! *************************************************************************
114 :
115 : !DEBUG
116 : !write(std_out,*)' MODULE predict_ga : enter '
117 : !ENDDEBUG
118 :
119 : ! gen dimension
120 :
121 4 : ndimen=3*natom
122 4 : nspinor=results_img(nimage,itimimage_eff)%nsppol
123 :
124 16 : ABI_MALLOC(coor,(ndimen,nimage))
125 12 : ABI_MALLOC(coor_old,(ndimen,nimage))
126 12 : ABI_MALLOC(acell,(3,nimage))
127 8 : ABI_MALLOC(acell_old,(3,nimage))
128 12 : ABI_MALLOC(rprim,(3,3,nimage))
129 8 : ABI_MALLOC(rprimd,(3,3,nimage))
130 8 : ABI_MALLOC(rprim_old,(3,3,nimage))
131 :
132 12 : ABI_MALLOC(zperm1,(natom))
133 8 : ABI_MALLOC(zperm2,(natom))
134 12 : ABI_MALLOC(ieperm,(nimage))
135 8 : ABI_MALLOC(ihperm,(nimage))
136 16 : ABI_MALLOC(ibgperm,(nspinor,nimage))
137 12 : ABI_MALLOC(ibgoptperm,(nspinor,nimage))
138 :
139 12 : ABI_MALLOC(etotal_img,(nimage))
140 8 : ABI_MALLOC(enthalpy_img,(nimage))
141 :
142 : !to store locally the calculated band gaps
143 16 : ABI_MALLOC(bg_img,(nspinor,nimage))
144 12 : ABI_MALLOC(bg_opt_img,(nspinor,nimage))
145 :
146 8 : ABI_MALLOC(fitness,(nimage))
147 :
148 12 : ABI_MALLOC(zcoor1,(natom))
149 8 : ABI_MALLOC(zcoor2,(natom))
150 12 : ABI_MALLOC(vson1,(ndimen))
151 8 : ABI_MALLOC(vson2,(ndimen))
152 :
153 44 : call initialize_perm(ieperm,nimage)
154 8 : call initialize_perm(ihperm,nimage)
155 8 : do ii=1,nspinor
156 8 : call initialize_perm(ibgperm(ii,:),nimage)
157 12 : call initialize_perm(ibgoptperm(ii,:),nimage)
158 : enddo
159 :
160 : !from gs_results_image take energies, rprim and acell and build energy and enthalpy vectors reordered from larger to smaller
161 :
162 44 : do iimage=1,nimage
163 40 : etotal_img(iimage)=results_img(iimage,itimimage_eff)%results_gs%etotal
164 40 : call convert_coortogen(results_img(iimage,itimimage_eff)%xred(:,:),coor_old(:,iimage),natom)
165 160 : acell_old(:,iimage)=results_img(iimage,itimimage_eff)%acell(:)
166 520 : rprim_old(:,:,iimage)=results_img(iimage,itimimage_eff)%rprim(:,:)
167 40 : if (results_img(iimage,itimimage_eff)%results_gs%gaps(3,1) > 0.0_dp) then
168 80 : bg_img(:,iimage)=results_img(iimage,itimimage_eff)%results_gs%gaps(1,:)
169 80 : bg_opt_img(:,iimage)=results_img(iimage,itimimage_eff)%results_gs%gaps(2,:)
170 : else
171 0 : bg_img(:,iimage)=-100.0_dp
172 0 : bg_opt_img(:,iimage)=-100.0_dp
173 : endif
174 40 : call mkrdim(acell_old(:,iimage),rprim_old(:,:,iimage),rprimd(:,:,iimage))
175 40 : call metric(gmet,gprimd,-1,rmet,rprimd(:,:,iimage),ucvol)
176 : enthalpy_img(iimage)=etotal_img(iimage) &
177 164 : & +sum(results_img(iimage,itimimage_eff)%results_gs%strten(1:3))*ucvol
178 : enddo
179 :
180 : !sort energies
181 :
182 4 : call sort_dp(nimage,etotal_img,ieperm,tol9)
183 :
184 : !sort enthalpies
185 :
186 4 : call sort_dp(nimage,enthalpy_img,ihperm,tol9)
187 :
188 : !sort band gaps
189 :
190 8 : do ii=1,nspinor
191 4 : call sort_dp(nimage,bg_img(ii,:),ibgperm(ii,:),tol9)
192 8 : call sort_dp(nimage,bg_opt_img(ii,:),ibgoptperm(ii,:),tol9)
193 : enddo
194 :
195 : ! Fitness is calculated
196 :
197 44 : sumH=0.d0
198 : Hmin=minval(enthalpy_img)
199 : Hmax=maxval(enthalpy_img)
200 44 : fitness = zero
201 :
202 4 : select case (ga_param%ga_fitness)
203 : case ( 1 )
204 :
205 : ! Function weighted over the difference with respect to the minimum
206 :
207 44 : do iimage=1,nimage
208 44 : sumh=sumh+enthalpy_img(iimage); fitness(iimage)=sumh
209 : enddo
210 4 : denom=nimage*enthalpy_img(nimage)-sumh
211 44 : do iimage=1,nimage
212 44 : fitness(iimage)=(iimage*enthalpy_img(nimage)-fitness(iimage))/denom
213 : enddo
214 :
215 : case ( 2 )
216 :
217 : ! Botzmann like, using the enthalpy difference.
218 :
219 0 : do iimage=1,nimage
220 0 : fitness(iimage)=exp(-one*(enthalpy_img(iimage)-enthalpy_img(1)))
221 0 : sumH = sumH + fitness(iimage)
222 : enddo
223 0 : fitness = fitness/sumH
224 0 : do iimage=2,nimage
225 0 : fitness(iimage)=fitness(iimage)+fitness(iimage-1)
226 : enddo
227 :
228 : case ( 3 )
229 :
230 : ! weighted ove the position in the ordered list, with probability 1/i (non uniform)
231 :
232 0 : do iimage=1,nimage
233 0 : fitness(iimage)=one/float(iimage)
234 0 : sumH = sumH + fitness(iimage)
235 : enddo
236 0 : fitness = fitness/sumH
237 4 : do iimage=2,nimage
238 0 : fitness(iimage)=fitness(iimage)+fitness(iimage-1)
239 : enddo
240 :
241 : end select
242 :
243 : ! do a single GA boocle
244 :
245 4 : indiv=0
246 :
247 : ! Selection over the best ga_opt_percent of the population
248 :
249 4 : nsurvivor=int(ga_param%ga_opt_percent*nimage)
250 :
251 4 : if (nsurvivor < one) nsurvivor=1
252 :
253 : ! pass coordinates,rprim, acell of survivors to next generation
254 :
255 12 : do iimage=1,nsurvivor
256 8 : indiv=indiv+1
257 200 : coor(:,iimage)=coor_old(:,ihperm(iimage))
258 32 : acell(:,iimage)=acell_old(:,ihperm(iimage))
259 108 : rprim(:,:,iimage)=rprim_old(:,:,ihperm(iimage))
260 : enddo
261 :
262 : ! complete the number of individuals of the generation by choosing them through GA
263 :
264 62 : do while(indiv<nimage)
265 :
266 : ! ga_n_rules corresponds to the number of chosen Genetic rules
267 :
268 58 : oper=ga_param%ga_rules(int(ga_param%ga_n_rules*uniformrandom(idum)+1))
269 :
270 10 : select case(oper)
271 : case(1) ! cut and splice elements after a randomization
272 10 : iparent1=choosefather(fitness,nimage,idum)
273 10 : iparent2=choosefather(fitness,nimage,idum)
274 10 : call convert_gentocoor(parent1,coor_old(:,ihperm(iparent1)),natom)
275 10 : call convert_gentocoor(parent2,coor_old(:,ihperm(iparent2)),natom)
276 : ! randomize the cell: random rotation and translation
277 10 : call randomize_parent(parent1,natom,idum)
278 10 : call randomize_parent(parent2,natom,idum)
279 : ! choose direction of cutting plane
280 10 : itmp = int(3*uniformrandom(idum)+1)
281 : ! order coordinates from small to large along that random direction axis of both parents
282 90 : zcoor1(:)=parent1(itmp,:)
283 90 : zcoor2(:)=parent2(itmp,:)
284 90 : call initialize_perm(zperm1,natom)
285 10 : call initialize_perm(zperm2,natom)
286 10 : call sort_dp(natom,zcoor1,zperm1,tol9)
287 10 : call sort_dp(natom,zcoor2,zperm2,tol9)
288 : ! choose the atom position to take the cut and take atoms below from one parent and above from the other parent
289 10 : itmp=int(natom*uniformrandom(idum)+1)
290 68 : do ii=1,itmp
291 232 : son1(:,ii)=parent1(:,zperm1(ii))
292 242 : son2(:,ii)=parent2(:,zperm2(ii))
293 : enddo
294 32 : do ii=itmp+1,natom
295 88 : son1(:,ii)=parent2(:,zperm2(ii))
296 98 : son2(:,ii)=parent1(:,zperm1(ii))
297 : enddo
298 : ! random combination of rprimd from parents
299 10 : call mkrdim(acell_old(:,ihperm(iparent1)),rprim_old(:,:,ihperm(iparent1)),rprimdparent1)
300 10 : call mkrdim(acell_old(:,ihperm(iparent2)),rprim_old(:,:,ihperm(iparent2)),rprimdparent2)
301 40 : do ii=1,3
302 130 : do jj=1,3
303 120 : rtmp(ii,jj)=uniformrandom(idum)
304 : enddo
305 : enddo
306 130 : rprimdson1=(one-rtmp)*rprimdparent1+rtmp*rprimdparent2
307 40 : do ii=1,3
308 130 : do jj=1,3
309 120 : rtmp(ii,jj)=uniformrandom(idum)
310 : enddo
311 : enddo
312 130 : rprimdson2=(one-rtmp)*rprimdparent1+rtmp*rprimdparent2
313 : !create acell and rprim from rprimd of springoffs
314 10 : call mkradim(acellson1,rprimson1,rprimdson1)
315 10 : call mkradim(acellson2,rprimson2,rprimdson2)
316 : ! check distances of atoms of every springoff
317 10 : if (checkatomicdist(natom,son1,rprimdson1)==0 .and. indiv<nimage) then
318 2 : indiv=indiv+1
319 : !if any fix coordinate restore the parent coordinate without modifications
320 18 : do ii=1,natom
321 66 : do jj=1,3
322 64 : if (ga_param%ga_iatfix(jj,ii) == 1) son1(jj,ii)=parent1(jj,ii)
323 : enddo
324 : enddo
325 2 : call convert_coortogen(son1,coor(:,indiv),natom)
326 8 : acell(:,indiv)=acellson1
327 26 : rprim(:,:,indiv)=rprimson1
328 : endif
329 10 : if (checkatomicdist(natom,son2,rprimdson2)==0 .and. indiv<nimage) then
330 6 : indiv=indiv+1
331 : !if any fix coordinate restore the parent coordinate without modifications
332 54 : do ii=1,natom
333 198 : do jj=1,3
334 192 : if (ga_param%ga_iatfix(jj,ii) == 1) son2(jj,ii)=parent1(jj,ii)
335 : enddo
336 : enddo
337 6 : call convert_coortogen(son2,coor(:,indiv),natom)
338 24 : acell(:,indiv)=acellson2
339 78 : rprim(:,:,indiv)=rprimson2
340 : endif
341 : case(2)! vector flip mutation
342 19 : iparent1=choosefather(fitness,nimage,idum)
343 19 : ii=int(ndimen*uniformrandom(idum)+1)
344 19 : jj=int(ndimen*uniformrandom(idum)+1)
345 19 : if (ii>jj) then
346 12 : call swap(ii,jj)
347 : end if
348 195 : vson1(1:ii)=coor_old(1:ii,ihperm(iparent1))
349 172 : if (jj<ndimen) vson1(jj+1:ndimen)=coor_old(jj+1:ndimen,ihperm(iparent1))
350 165 : do kk=ii,jj
351 165 : vson1(kk)=coor_old(jj+1-kk,ihperm(iparent1))
352 : enddo
353 247 : rprimson1=rprim_old(:,:,ihperm(iparent1))
354 76 : acellson1=acell_old(:,ihperm(iparent1))
355 19 : call convert_gentocoor(son1,vson1,natom)
356 19 : call mkrdim(acellson1,rprimson1,rprimdson1)
357 : !if any fix coordinate restore the parent coordinate without modifications
358 171 : do ii=1,natom
359 627 : do jj=1,3
360 608 : if (ga_param%ga_iatfix(jj,ii) == 1) son1(jj,ii)=parent1(jj,ii)
361 : enddo
362 : enddo
363 19 : call convert_coortogen(son1,vson1,natom)
364 19 : if (checkatomicdist(natom,son1,rprimdson1)==0 .and. indiv<nimage) then
365 7 : indiv=indiv+1
366 175 : coor(:,indiv)=vson1
367 28 : acell(:,indiv)=acellson1
368 91 : rprim(:,:,indiv)=rprimson1
369 : endif
370 : case(3) ! random strain - nondiagonal
371 12 : iparent1=choosefather(fitness,nimage,idum)
372 12 : rtmp(1,1)=one+gaussian_random(idum,0.1_dp)
373 12 : rtmp(2,2)=one+gaussian_random(idum,0.1_dp)
374 12 : rtmp(3,3)=one+gaussian_random(idum,0.1_dp)
375 12 : rtmp(1,2)=gaussian_random(idum,0.1_dp)*half
376 12 : rtmp(1,3)=gaussian_random(idum,0.1_dp)*half
377 12 : rtmp(2,3)=gaussian_random(idum,0.1_dp)*half
378 12 : rtmp(2,1)=rtmp(1,2)
379 12 : rtmp(3,1)=rtmp(1,3)
380 12 : rtmp(3,2)=rtmp(2,3)
381 480 : rprimson1=matmul(rtmp,rprim_old(:,:,ihperm(iparent1)))
382 312 : vson1=coor_old(:,ihperm(iparent1))
383 48 : acellson1=acell_old(:,ihperm(iparent1))
384 12 : call convert_gentocoor(son1,vson1,natom)
385 12 : call mkrdim(acellson1,rprimson1,rprimdson1)
386 : !if any fix coordinate restore the parent coordinate without modifications
387 108 : do ii=1,natom
388 396 : do jj=1,3
389 384 : if (ga_param%ga_iatfix(jj,ii) == 1) son1(jj,ii)=parent1(jj,ii)
390 : enddo
391 : enddo
392 12 : call convert_coortogen(son1,vson1,natom)
393 12 : if (checkatomicdist(natom,son1,rprimdson1)==0 .and. indiv<nimage) then
394 11 : indiv=indiv+1
395 275 : coor(:,indiv)=vson1
396 44 : acell(:,indiv)=acellson1
397 143 : rprim(:,:,indiv)=rprimson1
398 : endif
399 : case(4) ! coordinates mutation
400 17 : iparent1=choosefather(fitness,nimage,idum)
401 442 : vson1=coor_old(:,ihperm(iparent1))
402 17 : itmp=ndimen/4
403 17 : if (itmp<1) itmp=1
404 119 : do jj=1,itmp
405 102 : ii=int(ndimen*uniformrandom(idum)+1)
406 102 : vson1(ii)=vson1(ii)+0.15*uniformrandom(idum)
407 102 : if (vson1(ii)>one) vson1(ii)=vson1(ii)-one
408 119 : if (vson1(ii)<zero) vson1(ii)=vson1(ii)+one
409 : enddo
410 221 : rprimson1=rprim_old(:,:,ihperm(iparent1))
411 68 : acellson1=acell_old(:,ihperm(iparent1))
412 17 : call convert_gentocoor(son1,vson1,natom)
413 17 : call mkrdim(acellson1,rprimson1,rprimdson1)
414 : !if any fix coordinate restore the parent coordinate without modifications
415 153 : do ii=1,natom
416 561 : do jj=1,3
417 544 : if (ga_param%ga_iatfix(jj,ii) == 1) son1(jj,ii)=parent1(jj,ii)
418 : enddo
419 : enddo
420 17 : call convert_coortogen(son1,vson1,natom)
421 75 : if (checkatomicdist(natom,son1,rprimdson1)==0 .and. indiv<nimage) then
422 6 : indiv=indiv+1
423 150 : coor(:,indiv)=vson1
424 24 : acell(:,indiv)=acellson1
425 78 : rprim(:,:,indiv)=rprimson1
426 : endif
427 : end select
428 : enddo
429 :
430 4 : next_itimimage=itimimage_eff+1
431 4 : if (next_itimimage>ntimimage_stored) next_itimimage=1
432 :
433 44 : do iimage=1,nimage
434 200 : results_img(iimage,next_itimimage)%acell = acell(:,iimage)
435 560 : results_img(iimage,next_itimimage)%rprim = rprim(:,:,iimage)
436 1360 : results_img(iimage,next_itimimage)%vel = results_img(iimage,itimimage_eff)%vel
437 560 : results_img(iimage,next_itimimage)%vel_cell = results_img(iimage,itimimage_eff)%vel_cell
438 44 : call convert_gentocoor(results_img(iimage,next_itimimage)%xred,coor(:,iimage),natom)
439 : enddo
440 :
441 4 : ABI_FREE(coor)
442 4 : ABI_FREE(acell)
443 4 : ABI_FREE(acell_old)
444 4 : ABI_FREE(rprim)
445 4 : ABI_FREE(rprimd)
446 4 : ABI_FREE(rprim_old)
447 4 : ABI_FREE(zcoor1)
448 4 : ABI_FREE(zcoor2)
449 4 : ABI_FREE(coor_old)
450 4 : ABI_FREE(zperm1)
451 4 : ABI_FREE(zperm2)
452 4 : ABI_FREE(ieperm)
453 4 : ABI_FREE(ihperm)
454 4 : ABI_FREE(ibgperm)
455 4 : ABI_FREE(ibgoptperm)
456 4 : ABI_FREE(etotal_img)
457 4 : ABI_FREE(bg_img)
458 4 : ABI_FREE(bg_opt_img)
459 4 : ABI_FREE(enthalpy_img)
460 4 : ABI_FREE(fitness)
461 4 : ABI_FREE(vson1)
462 4 : ABI_FREE(vson2)
463 :
464 4 : end subroutine predict_ga
465 :
466 : !!*************** local subroutines
467 :
468 68 : INTEGER FUNCTION choosefather(fitf,n,idum)
469 :
470 : integer,intent(in) :: n
471 : integer,intent(inout) :: idum
472 : real(dp), dimension(:), intent(in) :: fitf
473 :
474 : real(dp) :: x1
475 : integer :: ii
476 :
477 68 : x1=uniformrandom(idum);
478 68 : choosefather=1
479 415 : do ii=2,n
480 415 : if(fitf(ii-1)<x1.and.x1<=fitf(ii))then
481 : choosefather=ii
482 : exit
483 : endif
484 : enddo
485 68 : end FUNCTION choosefather
486 :
487 : !!
488 :
489 12 : SUBROUTINE swap(a,b)
490 :
491 : integer, intent(inout) :: a,b
492 : integer :: dum
493 12 : dum=a; a=b; b=dum
494 12 : END SUBROUTINE swap
495 :
496 : !!
497 :
498 : INTEGER FUNCTION comp_indiv(distances,indiv,natom,nimage)
499 :
500 : !! comparing individuals from the same generation and check they
501 : !! are not too close. We compare all individuals with individual: indiv.
502 : !! We assume, all distances for a given individual are ordered and
503 : !! we define a metric from the total difference distances between individuals.
504 : !! if comp_indiv is 0, means that two individuals are two close.
505 :
506 : integer, intent(in) :: indiv,natom,nimage
507 : real(dp),intent(in) ::distances(natom,nimage)
508 : real(dp) :: diff
509 :
510 : integer :: ii,jj
511 :
512 : comp_indiv=1
513 :
514 : do ii=1,indiv-1
515 : diff=0.0_dp
516 : do jj=1,natom
517 : diff=(distances(jj,indiv)-distances(jj,ii))**2
518 : enddo
519 : if (diff.le.0.001_dp) comp_indiv=0
520 : enddo
521 :
522 : end FUNCTION comp_indiv
523 :
524 : !!
525 :
526 20 : SUBROUTINE randomize_parent(parent,natom,idum)
527 :
528 : ! take a parent and randomize the positions of the atoms
529 :
530 : integer,intent(in) :: natom
531 : integer,intent(inout) :: idum
532 : real(dp),intent(inout) :: parent(3,natom)
533 :
534 : real(dp) :: tmp(3),rot(3,3)
535 : integer :: ii,jj
536 :
537 : !random rotation
538 :
539 180 : do ii=1,natom
540 160 : if (uniformrandom(idum)>half) then
541 312 : do jj=1,3
542 312 : tmp(jj)=uniformrandom(idum)*two_pi
543 : enddo
544 78 : rot(1,1)=cos(tmp(3))*cos(tmp(2))
545 78 : rot(1,2)=cos(tmp(3))*sin(tmp(1))*sin(tmp(2))-cos(tmp(1))*sin(tmp(3))
546 78 : rot(1,3)=cos(tmp(1))*cos(tmp(3))*sin(tmp(2))+sin(tmp(1))*sin(tmp(3))
547 78 : rot(2,1)=cos(tmp(2))*sin(tmp(3))
548 78 : rot(2,2)=cos(tmp(1))*cos(tmp(3))+sin(tmp(1))*sin(tmp(2))*sin(tmp(3))
549 78 : rot(2,3)=cos(tmp(1))*sin(tmp(2))*sin(tmp(3))-cos(tmp(3))*sin(tmp(1))
550 78 : rot(3,1)=-sin(tmp(2))
551 78 : rot(3,2)=cos(tmp(2))*sin(tmp(1))
552 78 : rot(3,3)=cos(tmp(1))*cos(tmp(2))
553 1482 : parent(:,ii)=matmul(rot(:,:),parent(:,ii))
554 : endif
555 :
556 : ! random reflection
557 160 : if (uniformrandom(idum)>half) then
558 336 : parent(:,ii)=-parent(:,ii)
559 : endif
560 :
561 : ! fold back all coordinates to reduced coordinates in cube [0,1]
562 660 : parent(:,ii)=parent(:,ii)-anint((parent(:,ii)-one)/two)
563 :
564 : enddo
565 :
566 : ! shift all atoms along a random direction
567 :
568 80 : do ii=1,3
569 80 : tmp(ii)=uniformrandom(idum)
570 : enddo
571 :
572 180 : do ii=1,natom
573 640 : parent(:,ii)=parent(:,ii)+tmp(:)
574 660 : parent(:,ii)=parent(:,ii)-anint(parent(:,ii)/two)
575 : enddo
576 :
577 20 : end SUBROUTINE randomize_parent
578 :
579 : !!
580 :
581 108 : SUBROUTINE convert_gentocoor(parent,coor,natom)
582 :
583 : ! convert gene (single vector) to coordinates (3,natom)
584 :
585 : integer,intent(in) :: natom
586 : double precision,intent(out) :: parent(3,natom)
587 : double precision,intent(in) :: coor(3*natom)
588 :
589 : integer :: ii,jj
590 :
591 :
592 972 : do ii=1,natom
593 864 : jj=(ii-1)*3
594 3564 : parent(:,ii)=coor(jj+1:jj+3)
595 : enddo
596 :
597 108 : end SUBROUTINE convert_gentocoor
598 :
599 : !!
600 :
601 96 : SUBROUTINE convert_coortogen(parent,coor,natom)
602 :
603 : ! convert coordinates in gene notation (single vector with all coordinates)
604 :
605 : integer,intent(in) :: natom
606 : double precision,intent(in) :: parent(3,natom)
607 : double precision,intent(out) :: coor(3*natom)
608 :
609 : integer :: ii,jj
610 :
611 864 : do ii=1,natom
612 768 : jj=(ii-1)*3
613 3168 : coor(jj+1:jj+3)=parent(:,ii)
614 : enddo
615 :
616 96 : end SUBROUTINE convert_coortogen
617 :
618 : !!
619 :
620 8 : SUBROUTINE initialize_perm(iperm,nimage)
621 :
622 : !! initialize the vector iperm with corresponding indices.
623 :
624 : integer, intent(in) :: nimage
625 : integer, intent(inout) :: iperm(nimage)
626 : integer :: ii
627 :
628 :
629 356 : do ii=1,nimage
630 356 : iperm(ii)=ii
631 : enddo
632 :
633 : end SUBROUTINE initialize_perm
634 :
635 : ! if after a genetic rule, check if two atoms in the same gen
636 : ! are too close
637 :
638 :
639 68 : INTEGER FUNCTION checkatomicdist(natom,coord,rprimd)
640 :
641 : !! check if two atoms are two close.
642 : !! if they are, checkatomicdist=0
643 :
644 : integer, intent(in) :: natom
645 : real(dp), intent(in) :: coord(3,natom),rprimd(3,3)
646 : real(dp) :: d,v1(3),v2(3)
647 : integer :: ii,jj
648 :
649 68 : checkatomicdist=0
650 612 : do ii=1,natom
651 2176 : v1=coord(:,ii)
652 2391 : do jj=ii+1,natom
653 7288 : v2=coord(:,jj)
654 1822 : d=dist2(v1,v2,rprimd,1)
655 2323 : if (d<1.0d0) then
656 : checkatomicdist=1
657 : EXIT
658 : endif
659 : enddo
660 : enddo
661 :
662 68 : END FUNCTION checkatomicdist
663 :
664 72 : DOUBLE PRECISION FUNCTION gaussian_random(idum,sigma)
665 :
666 :
667 : integer,intent(inout) :: idum
668 : real(dp), intent(in) :: sigma
669 :
670 : real(dp) :: r1,r2,w
671 : !! The polar form of the Box-Muller transformation
672 72 : w=two
673 159 : do while (w >= one .or. w == zero)
674 87 : r1=two*uniformrandom(idum)-one
675 87 : r2=two*uniformrandom(idum)-one
676 87 : w=r1*r1+r2*r2
677 : enddo
678 72 : w=sqrt(-two*log(w)/w)
679 72 : gaussian_random=r2*w*sigma
680 :
681 72 : END FUNCTION gaussian_random
682 :
683 : end MODULE m_use_ga
|
