Branch data 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
|
