atomes 1.1.16
atomes: an atomic scale modeling tool box
Loading...
Searching...
No Matches
rings-king.F90
Go to the documentation of this file.
1! This file is part of the 'atomes' software.
2!
3! 'atomes' is free software: you can redistribute it and/or modify it under the terms
4! of the GNU Affero General Public License as published by the Free Software Foundation,
5! either version 3 of the License, or (at your option) any later version.
6!
7! 'atomes' is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
8! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
9! See the GNU General Public License for more details.
10!
11! You should have received a copy of the GNU Affero General Public License along with 'atomes'.
12! If not, see <https://www.gnu.org/licenses/>
13!
14! Copyright (C) 2022-2024 by CNRS and University of Strasbourg
15!
20
21SUBROUTINE setup_cpat_vpat_ring (NAT, STR, CONT, VOIS, CPT, VPT)
22
23USE parameters
24
25IMPLICIT NONE
26
27INTEGER, INTENT(IN) :: NAT, STR
28INTEGER, DIMENSION(NAT,NS), INTENT(IN):: CONT
29INTEGER, DIMENSION(MAXN,NAT,NS), INTENT(IN) :: VOIS
30INTEGER, DIMENSION(NAT), INTENT(INOUT):: CPT
31INTEGER, DIMENSION(NAT,MAXN), INTENT(INOUT) :: VPT
32INTEGER :: RAB, RAC, RAD
33
34!do RAB=1, NAT
35! write (6, '("At= ",i4," Neigh= ",i2)') RAB, CONTJ(RAB,STR)
36! do RAC=1, CONTJ(RAB,STR)
37! write (6, '(" i= ",i2," Vois= ",i4)') RAC, VOISJ(RAC,RAB,STR)
38! enddo
39!enddo
40
41cpt(:) = 0
42vpt(:,:) = 0
43
44do rab=1, nat
45 if (cont(rab,str) .gt. 1) then
46 rac = 0
47 do rad=1, cont(rab,str)
48 if (cont(vois(rad,rab,str),str) .gt. 1) then
49 rac=rac+1
50 vpt(rab,rac) = vois(rad,rab,str)
51 endif
52 enddo
53 if (rac .ge. 2) then
54 cpt(rab) = rac
55 endif
56 endif
57enddo
58
59END SUBROUTINE
60
61SUBROUTINE setup_cpat_vpat (STR, NAT)
62
63USE parameters
64
65IMPLICIT NONE
66
67INTEGER, INTENT(IN) :: STR, NAT
68INTEGER :: RAB, RAC, RAD
69
70!do RAB=1, NAT
71! write (6, '("At= ",i4," Neigh= ",i2)') RAB, CONTJ(RAB,STR)
72! do RAC=1, CONTJ(RAB,STR)
73! write (6, '(" i= ",i2," Vois= ",i4)') RAC, VOISJ(RAC,RAB,STR)
74! enddo
75!enddo
76
77
78do rab=1, nat
79 cpat(rab) = 0
80 if (contj(rab,str) .gt. 1) then
81 rac = 0
82 do rad=1, contj(rab,str)
83 if (contj(voisj(rad,rab,str),str) .gt. 1) then
84 rac=rac+1
85 vpat(rab,rac) = voisj(rad,rab,str)
86 endif
87 enddo
88 if (rac .ge. 2) then
89 cpat(rab) = rac
90 endif
91 endif
92enddo
93
94END SUBROUTINE
95
96INTEGER FUNCTION king_rings()
97
98USE parameters
99
100#ifdef OPENMP
101!$ USE OMP_LIB
102#endif
103IMPLICIT NONE
104
105INTEGER :: ar
106#ifdef OPENMP
107INTEGER :: numth
108LOGICAL :: doatoms
109#endif
110
111INTERFACE
112 INTEGER FUNCTION recrings(VID)
113 INTEGER, INTENT(IN) :: vid
114 END FUNCTION
115END INTERFACE
116
117king_rings = 0
118ar = 0
119if (.not.allrings) ar=1
120
121#ifdef OPENMP
122numth = omp_get_max_threads()
123doatoms=.false.
124! Dynamic allocation of pointers and tables used in the "RINGS" subroutine
125if (ns.ge.1 .and. ns.lt.numth) then
126 if (numth .ge. 2*(ns-1)) then
127 doatoms=.true.
128 else
129 numth=ns
130 endif
131endif
132
133if (all_atoms) doatoms=.true.
134
135if (doatoms) then
136 if (na.lt.numth) numth=na
137#ifdef DEBUG
138 write (6, *) "OpenMP on atoms, NUMTH= ",numth
139#endif
140 call king_ring_search_atoms (ar, numth)
141else
142#ifdef DEBUG
143 write (6, *) "OpenMP on MD steps, NUMTH= ",numth
144#endif
145 call king_ring_search_steps (ar, numth)
146endif
147#else
148call king_ring_search_steps (ar)
149#endif
150
151king_rings = recrings(ar)
152
153END FUNCTION
154
155#ifdef OPENMP
156SUBROUTINE king_ring_search_atoms (ARI, NUMTH)
157
158USE parameters
159
160!$ USE OMP_LIB
161IMPLICIT NONE
162
163INTEGER, INTENT(IN) :: ari, numth
164TYPE (ring), DIMENSION(:), ALLOCATABLE :: the_ring
165INTEGER, DIMENSION(:), ALLOCATABLE :: tring, indte, indth
166INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: savring, ordring
167INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: savr, ordr
168INTEGER :: lora, lorb, lorc, ri
169
170INTERFACE
171 RECURSIVE SUBROUTINE inside_ring (THE_RING, FND, S_IR, AI_IR, RID, TAE, TAH, LRA, LRB, &
172 NRPAT, RSAVED, OSAVED, TRING, INDE, INDH, RESL, CPT, VPT)
173 USE parameters
174 TYPE (ring), DIMENSION(TAILLD), INTENT(INOUT) :: the_ring
175 LOGICAL, INTENT(INOUT) :: fnd
176 INTEGER, INTENT(IN) :: s_ir, ai_ir, rid, lra, lrb
177 INTEGER, INTENT(INOUT) :: tae, tah
178 INTEGER, DIMENSION(NA), INTENT(INOUT) :: nrpat
179 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: rsaved, osaved
180 INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: inde, indh
181 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: resl, tring
182 INTEGER, DIMENSION(NA,MAXN), INTENT(IN) :: vpt
183 INTEGER, DIMENSION(NA), INTENT(IN):: cpt
184 END SUBROUTINE
185 INTEGER FUNCTION rings_to_ogl_menu (IDSEARCH, NRI)
186 USE parameters
187 INTEGER, INTENT(IN) :: idsearch
188 INTEGER, DIMENSION(TAILLR, NS), INTENT(IN) :: nri
189 END FUNCTION
190 INTEGER FUNCTION rings_to_ogl (STEP, IDSEARCH, NRI, RSAVED, OSAVED)
191 USE parameters
192 INTEGER, INTENT(IN) :: step, idsearch
193 INTEGER, DIMENSION(TAILLR,NS), INTENT(IN) :: nri
194 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(IN) :: rsaved, osaved
195 END FUNCTION
196 SUBROUTINE del_this_ring (TLED, RSAVED, OSAVED, TRING, RESL, INDT)
197 USE parameters
198 INTEGER, INTENT(IN) :: tled
199 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: rsaved, osaved
200 INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: indt
201 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: resl, tring
202 END SUBROUTINE
203END INTERFACE
204
205ri = 0
206if(allocated(savring)) deallocate(savring)
207allocate(savring(taillr,numa,taillr), stat=err)
208if (err .ne. 0) then
209 alc_tab="SAVRING"
210 alc=.true.
211 goto 001
212endif
213if(allocated(ordring)) deallocate(ordring)
214allocate(ordring(taillr,numa,taillr), stat=err)
215if (err .ne. 0) then
216 alc_tab="ORDRING"
217 alc=.true.
218 goto 001
219endif
220if(allocated(cpat)) deallocate(cpat)
221allocate(cpat(na), stat=err)
222if (err .ne. 0) then
223 alc_tab="CPAT"
224 alc=.true.
225 goto 001
226endif
227 if(allocated(vpat)) deallocate(vpat)
228allocate(vpat(na,maxn), stat=err)
229if (err .ne. 0) then
230 alc_tab="VPAT"
231 alc=.true.
232 goto 001
233endif
234
235do i=1, ns
236
237 p=0
238
239 savring(:,:,:)=0
240 ordring(:,:,:)=0
241 call setup_cpat_vpat_ring (na, i, contj, voisj, cpat, vpat)
242
243 ! OpenMP on atoms only
244 !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
245 !$OMP& PRIVATE(MAXAT, MINAT, RPAT, APNA, SAUT, RUNSEARCH, &
246 !$OMP& j, k, l, m, n, o, LORA, LORB, LORC, TAILLE, TAILLH, THE_RING, RES_LIST, &
247 !$OMP& FOUND, ERR, SAVR, ORDR, TRING, INDTE, INDTH) &
248 !$OMP& SHARED(i, p, NUMTH, NS, NA, TLT, NSP, LOT, TAILLR, TAILLD, CONTJ, VOISJ, CPAT, VPAT, &
249 !$OMP& NUMA, FACTATRING, ATRING, MAXPNA, MINPNA, DOAMPAT, AMPAT, ABAB, NO_HOMO, ALLRINGS, &
250 !$OMP& TBR, ALC, ALC_TAB, NCELLS, PBC, MAXN, SAVRING, ORDRING, NRING, INDRING, PNA, ri)
251
252 if (allocated(rpat)) deallocate(rpat)
253 allocate(rpat(na), stat=err)
254 if (err .ne. 0) then
255 alc_tab="RPAT"
256 alc=.true.
257 goto 002
258 endif
259 if(allocated(res_list)) deallocate(res_list)
260 allocate(res_list(taillr), stat=err)
261 if (err .ne. 0) then
262 alc_tab="RES_LIST"
263 alc=.true.
264 goto 002
265 endif
266 if(allocated(indte)) deallocate(indte)
267 allocate(indte(numa), stat=err)
268 if (err .ne. 0) then
269 alc_tab="INDTE"
270 goto 002
271 endif
272 if(allocated(indth)) deallocate(indth)
273 allocate(indth(numa), stat=err)
274 if (err .ne. 0) then
275 alc_tab="INDTH"
276 goto 002
277 endif
278 if(allocated(apna)) deallocate(apna)
279 allocate(apna(taillr), stat=err)
280 if (err .ne. 0) then
281 alc_tab="APNA"
282 alc=.true.
283 goto 002
284 endif
285 if(allocated(savr)) deallocate(savr)
286 allocate(savr(taillr,numa,taillr), stat=err)
287 if (err .ne. 0) then
288 alc_tab="SAVR"
289 alc=.true.
290 goto 002
291 endif
292 if(allocated(ordr)) deallocate(ordr)
293 allocate(ordr(taillr,numa,taillr), stat=err)
294 if (err .ne. 0) then
295 alc_tab="ORDR"
296 alc=.true.
297 goto 002
298 endif
299 if(allocated(tring)) deallocate(tring)
300 allocate(tring(taillr), stat=err)
301 if (err .ne. 0) then
302 alc_tab="TRING"
303 alc=.true.
304 goto 002
305 endif
306 if(allocated(the_ring)) deallocate(the_ring)
307 allocate(the_ring(tailld), stat=err)
308 if (err .ne. 0) then
309 alc_tab="THE_RING"
310 alc=.true.
311 goto 002
312 endif
313
314 tring(:)=0
315 savr(:,:,:)=0
316 ordr(:,:,:)=0
317 !$OMP DO SCHEDULE(STATIC,NA/NUMTH)
318 do j=1, na
319
320 if (tbr .or. alc) goto 003
321 if (tlt .eq. nsp+1 .or. lot(j) .eq. tlt) then
322
323 !$OMP CRITICAL
324 p = p+1
325 o = p
326 !$OMP END CRITICAL
327 apna(:) = 0
328 minat=tailld
329 maxat=1
330 saut=.true.
331
332 if (cpat(j).ge.2) then
333
334 do l=1, cpat(j)-1
335
336 do m=l+1, cpat(j)
337
338 lora=lot(j)
339 lorb=lot(vpat(j,l))
340 lorc=lot(vpat(j,m))
341 if (abab) then
342 if ((lorb .ne. lora) .and. (lorc .ne. lora) .and. (lorb .eq. lorc)) then
343 runsearch=.true.
344 else
345 runsearch=.false.
346 endif
347 else
348 runsearch=.true.
349 endif
350
351 if (runsearch) then
352
353 found=.false.
354 rpat(:) = 0
355 do k=1, cpat(j)
356 rpat(vpat(j,k)) = 1
357 enddo
358 rpat(vpat(j,l)) = 0
359 the_ring(1)%ATOM=vpat(j,l)
360 the_ring(1)%SPEC=lorb
361 the_ring(1)%NEIGHBOR=1
362 the_ring(2)%ATOM=j
363 the_ring(2)%SPEC=lora
364 the_ring(2)%NEIGHBOR=1
365 rpat(j)=1
366 taille=taillr
367 taillh=taillr
368 res_list(:) = 0
369 indte(:) = 0
370 indth(:) = 0
371 the_ring(3)%ATOM=vpat(j,m)
372 the_ring(3)%SPEC=lorc
373 the_ring(3)%NEIGHBOR=cpat(vpat(j,m))
374
375 call inside_ring (the_ring, found, i, vpat(j,m), 3, taille, taillh, lora, lorb, &
376 rpat, savr, ordr, tring, indte, indth, res_list, cpat, vpat)
377 if (alc) alc_tab="INSIDE_RING"
378 if (tbr .or. alc) goto 003
379
380 if (allrings) then
381 do k=3, taillr
382 do n=1, tring(k)
383 if (indte(n).ne.0 .and. apna(k).eq.0) then
384 apna(k)=1
385 minat=min(minat,k)
386 maxat=max(maxat,k)
387 ! if (FACTATPNA) ATPNA(k,o,i)=1
388 saut=.false.
389 endif
390 enddo
391 enddo
392 else if (found) then
393 if (apna(taille) .eq. 0) then
394 apna(taille)=1
395 minat=min(minat,taille)
396 maxat=max(maxat,taille)
397 ! if (FACTATPNA) ATPNA(TAILLE,o,i)=1
398 saut=.false.
399 endif
400 else if (no_homo .and. res_list(taillh) .ne. 0) then
401 ! Some shortest ring with HP bonds and < TAILLR were found
402 ! delete these rings now
403 call del_this_ring (taillh, savr, ordr, tring, res_list, indth)
404 else if (doampat) then
405 if (contj(vpat(j,l),i) .ge. 2 .and. contj(vpat(j,m),i) .ge. 2) then
406 !$OMP ATOMIC
407 ampat(o,i)=ampat(o,i)+1
408 endif
409 endif
410 endif
411 enddo
412 enddo
413
414 endif
415
416 if (.not. saut) then
417 do k=3, taillr
418 do l=3, taillr
419 if (apna(k).eq.1 .and. apna(l).eq.1) then
420 !$OMP ATOMIC
421 pna(k,l,i)=pna(k,l,i)+1
422 endif
423 enddo
424 enddo
425 !$OMP ATOMIC
426 maxpna(maxat,i)=maxpna(maxat,i)+1
427 !$OMP ATOMIC
428 minpna(minat,i)=minpna(minat,i)+1
429 endif
430 endif
431
432 003 continue
433 enddo
434 !$OMP END DO NOWAIT
435 if (tbr .or. alc) goto 002
436 !$OMP CRITICAL
437 do k=3, taillr
438 if (tring(k).gt.0) then
439 if (nring(k,i).gt.0) then
440 o = 0
441 do l=1, tring(k)
442 do m=1, nring(k,i)
443 saut=.true.
444 do n=1, k
445 if (savring(k,m,n) .ne. savr(k,l,n)) then
446 saut=.false.
447 exit
448 endif
449 enddo
450 if (saut) exit
451 enddo
452 if (.not.saut) then
453 o = o + 1
454 if (nring(k,i)+o .gt. numa) then
455 tbr=.true.
456 goto 004
457 endif
458 do m=1, k
459 savring(k,nring(k,i)+o,m) = savr(k,l,m)
460 ordring(k,nring(k,i)+o,m) = ordr(k,l,m)
461 enddo
462 endif
463 enddo
464 nring(k,i)=nring(k,i)+o
465 else
466 do l=1, tring(k)
467 do m=1, k
468 savring(k,l,m) = savr(k,l,m)
469 ordring(k,l,m) = ordr(k,l,m)
470 enddo
471 enddo
472 nring(k,i) = tring(k)
473 endif
474 endif
475 enddo
476
477 004 continue
478 !$OMP END CRITICAL
479
480 002 continue
481
482 if (allocated(rpat)) deallocate (rpat)
483 if (allocated(res_list)) deallocate (res_list)
484 if (allocated(indte)) deallocate (indte)
485 if (allocated(indth)) deallocate (indth)
486 if (allocated(apna)) deallocate (apna)
487 if (allocated(tring)) deallocate (tring)
488 if (allocated(savr)) deallocate (savr)
489 if (allocated(ordr)) deallocate (ordr)
490 if (allocated(the_ring)) deallocate (the_ring)
491 !$OMP END PARALLEL
492 if (alc .or. tbr) goto 001
493
494 !do j=3, TAILLR
495 ! write (6, '("s= ",i4,", j= ",i2,", nr(",i2,",",i4,")= ",i2)') i,j,j,i, NRING(j,i)
496 ! if (NRING(j,i) .gt. 0) then
497 ! do k=1, NRING(j,i)
498 ! write (6, *) " k= ",k,", R(o)= ",ORDRING(j,k,1:j)
499 ! enddo
500 ! endif
501 !enddo
502 ri = ri + rings_to_ogl(i, ari, nring, savring, ordring)
503
504enddo
505
506001 continue
507
508if (alc) then
509 call show_error ("Impossible to allocate memory"//char(0), &
510 "Subroutine: KING_RING_SEARCH_ATOMS"//char(0), "Table: "//alc_tab(1:len_trim(alc_tab))//char(0))
511endif
512
513if (allocated(cpat)) deallocate (cpat)
514if (allocated(vpat)) deallocate (vpat)
515if (allocated(savring)) deallocate (savring)
516if (allocated(ordring)) deallocate (ordring)
517
518if (ri .eq. ns) ri = rings_to_ogl_menu(ari, nring)
519
520END SUBROUTINE
521#endif
522
523#ifdef OPENMP
524SUBROUTINE king_ring_search_steps (ARI, NUMTH)
525
526USE parameters
527!$ USE OMP_LIB
528IMPLICIT NONE
529INTEGER, INTENT(IN) :: ari, numth
530#else
531SUBROUTINE king_ring_search_steps (ARI)
532
533USE parameters
534IMPLICIT NONE
535INTEGER, INTENT(IN) :: ARI
536#endif
537TYPE (RING), DIMENSION(:), ALLOCATABLE :: THE_RING
538INTEGER, DIMENSION(:), ALLOCATABLE :: TRING, INDTE, INDTH
539INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: SAVRING, ORDRING
540INTEGER :: LORA, LORB, LORC, ri
541
542INTERFACE
543 RECURSIVE SUBROUTINE inside_ring (THE_RING, FND, S_IR, AI_IR, RID, TAE, TAH, LRA, LRB, &
544 NRPAT, RSAVED, OSAVED, TRING, INDE, INDH, RESL, CPT, VPT)
545 USE parameters
546 TYPE (RING), DIMENSION(TAILLD), INTENT(INOUT) :: THE_RING
547 LOGICAL, INTENT(INOUT) :: FND
548 INTEGER, INTENT(IN) :: S_IR, AI_IR, RID
549 INTEGER, INTENT(INOUT) :: TAE, TAH
550 INTEGER, INTENT(IN) :: LRA, LRB
551 INTEGER, DIMENSION(NA), INTENT(INOUT) :: NRPAT
552 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: RSAVED, OSAVED
553 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: TRING
554 INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: INDE, INDH
555 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: RESL
556 INTEGER, DIMENSION(NA), INTENT(IN):: CPT
557 INTEGER, DIMENSION(NA,MAXN), INTENT(IN) :: VPT
558 END SUBROUTINE
559 INTEGER FUNCTION rings_to_ogl_menu (IDSEARCH, NRI)
560 USE parameters
561 INTEGER, INTENT(IN) :: IDSEARCH
562 INTEGER, DIMENSION(TAILLR, NS), INTENT(IN) :: NRI
563 END FUNCTION
564 INTEGER FUNCTION rings_to_ogl (STEP, IDSEARCH, NRI, RSAVED, OSAVED)
565 USE parameters
566 INTEGER, INTENT(IN) :: STEP, IDSEARCH
567 INTEGER, DIMENSION(TAILLR,NS), INTENT(IN) :: NRI
568 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(IN) :: RSAVED, OSAVED
569 END FUNCTION
570 SUBROUTINE del_this_ring (TLED, RSAVED, OSAVED, TRING, RESL, INDT)
571 USE parameters
572 INTEGER, INTENT(IN) :: TLED
573 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: RSAVED, OSAVED
574 INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: INDT
575 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: RESL, TRING
576 END SUBROUTINE
577END INTERFACE
578
579ri = 0
580#ifdef OPENMP
581! OpenMP on steps only
582!$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
583!$OMP& PRIVATE(TAILLE, TAILLH, MAXAT, MINAT, SAUT, RUNSEARCH, &
584!$OMP& i, j, k, l, m, n, o, p, LORA, LORB, LORC, THE_RING, RES_LIST, INDTE, INDTH, APNA, &
585!$OMP& FOUND, ERR, TRING, SAVRING, ORDRING, RPAT, CPAT, VPAT) &
586!$OMP& SHARED(ARI, NUMTH, NS, NA, TLT, NSP, LOT, TAILLR, TAILLD, CONTJ, VOISJ, &
587!$OMP& NUMA, FACTATRING, ATRING, MAXPNA, MINPNA, DOAMPAT, AMPAT, ABAB, NO_HOMO, ALLRINGS, &
588!$OMP& TBR, ALC, ALC_TAB, NCELLS, THE_BOX, FULLPOS, PBC, MAXN, NRING, INDRING, PNA, ri)
589#endif
590
591if(allocated(savring)) deallocate(savring)
592allocate(savring(taillr,numa,taillr), stat=err)
593if (err .ne. 0) then
594 alc_tab="SAVRING"
595 alc=.true.
596 goto 001
597endif
598if(allocated(ordring)) deallocate(ordring)
599allocate(ordring(taillr,numa,taillr), stat=err)
600if (err .ne. 0) then
601 alc_tab="ORDRING"
602 alc=.true.
603 goto 001
604endif
605if(allocated(cpat)) deallocate(cpat)
606allocate(cpat(na), stat=err)
607if (err .ne. 0) then
608 alc_tab="CPAT"
609 alc=.true.
610 goto 001
611endif
612if(allocated(vpat)) deallocate(vpat)
613allocate(vpat(na,maxn), stat=err)
614if (err .ne. 0) then
615 alc_tab="VPAT"
616 alc=.true.
617 goto 001
618endif
619if (allocated(rpat)) deallocate(rpat)
620allocate(rpat(na), stat=err)
621if (err .ne. 0) then
622 alc_tab="RPAT"
623 alc=.true.
624 goto 001
625endif
626if(allocated(res_list)) deallocate(res_list)
627allocate(res_list(taillr), stat=err)
628if (err .ne. 0) then
629 alc_tab="RES_LIST"
630 alc=.true.
631 goto 001
632endif
633if(allocated(indte)) deallocate(indte)
634allocate(indte(numa), stat=err)
635if (err .ne. 0) then
636 alc_tab="INDTE"
637 goto 001
638endif
639if(allocated(indth)) deallocate(indth)
640allocate(indth(numa), stat=err)
641if (err .ne. 0) then
642 alc_tab="INDTH"
643 goto 001
644endif
645if(allocated(apna)) deallocate(apna)
646allocate(apna(taillr), stat=err)
647if (err .ne. 0) then
648 alc_tab="APNA"
649 alc=.true.
650 goto 001
651endif
652if(allocated(tring)) deallocate(tring)
653allocate(tring(taillr), stat=err)
654if (err .ne. 0) then
655 alc_tab="TRING"
656 alc=.true.
657 goto 001
658endif
659if(allocated(the_ring)) deallocate(the_ring)
660allocate(the_ring(tailld), stat=err)
661if (err .ne. 0) then
662 alc_tab="THE_RING"
663 alc=.true.
664 goto 001
665endif
666
667#ifdef OPENMP
668!$OMP DO SCHEDULE(STATIC,NS/NUMTH)
669#endif
670do i=1, ns
671
672 if (tbr .or. alc) goto 002
673 savring(:,:,:)=0
674 ordring(:,:,:)=0
675 tring(:)=0
676 call setup_cpat_vpat_ring (na, i, contj, voisj, cpat, vpat)
677
678 o=0
679 do j=1, na
680
681 if (tlt .eq. nsp+1 .or. lot(j) .eq. tlt) then
682
683 o=o+1
684 apna(:) = 0
685 minat=taillr
686 maxat=1
687 saut=.true.
688
689 if (cpat(j).ge.2) then
690
691 do l=1, cpat(j)-1
692
693 do m=l+1, cpat(j)
694
695 lora=lot(j)
696 lorb=lot(vpat(j,l))
697 lorc=lot(vpat(j,m))
698 if (abab) then
699 if ((lorb .ne. lora) .and. (lorc .ne. lora) .and. (lorb .eq. lorc)) then
700 runsearch=.true.
701 else
702 runsearch=.false.
703 endif
704 else
705 runsearch=.true.
706 endif
707
708 if (runsearch) then
709
710 found=.false.
711 rpat(:) = 0
712 do k=1, cpat(j)
713 rpat(vpat(j,k)) = 1
714 enddo
715 p = vpat(j,l)
716 rpat(p) = 0
717 rpat(j)=1
718 the_ring(1)%ATOM=p
719 the_ring(1)%SPEC=lorb
720 the_ring(1)%NEIGHBOR=1
721 the_ring(2)%ATOM=j
722 the_ring(2)%SPEC=lora
723 the_ring(2)%NEIGHBOR=1
724 the_ring(3)%ATOM=vpat(j,m)
725 the_ring(3)%SPEC=lorc
726 the_ring(3)%NEIGHBOR=cpat(vpat(j,m))
727 taille=taillr
728 taillh=taillr
729 res_list(:) = 0
730 indte(:) = 0
731 indth(:) = 0
732 call inside_ring (the_ring, found, i, vpat(j,m), 3, taille, taillh, lora, lorb, &
733 rpat, savring, ordring, tring, indte, indth, res_list, cpat, vpat)
734 if (alc) alc_tab="INSIDE_RING"
735 if (tbr .or. alc) goto 002
736
737 if (allrings) then
738 do k=3, taillr
739 do n=1, tring(k)
740 if (indte(n).ne.0 .and. apna(k).eq.0) then
741 apna(k)=1
742 minat=min(minat,k)
743 maxat=max(maxat,k)
744 ! if (FACTATPNA) ATPNA(k,o,i)=1
745 saut=.false.
746 endif
747 enddo
748 enddo
749 else if (found) then
750 if (apna(taille) .eq. 0) then
751 apna(taille)=1
752 minat=min(minat,taille)
753 maxat=max(maxat,taille)
754 ! if (FACTATPNA) ATPNA(TAILLE,o,i)=1
755 saut=.false.
756 endif
757 else if (no_homo .and. res_list(taillh) .ne. 0) then
758 ! Some shortest ring with HP bonds and < TAILLR were found
759 ! delete these rings now
760 call del_this_ring (taillh, savring, ordring, tring, res_list, indth)
761 else if (doampat) then
762 if (contj(vpat(j,l),i) .ge. 2 .and. contj(vpat(j,m),i) .ge. 2) ampat(o,i)=ampat(o,i)+1
763 endif
764 endif
765 enddo
766 enddo
767
768 endif
769
770 if (.not. saut) then
771 do k=3, taillr
772 do l=3, taillr
773 if (apna(k).eq.1 .and. apna(l).eq.1) then
774 pna(k,l,i)=pna(k,l,i)+1
775 endif
776 enddo
777 enddo
778 maxpna(maxat,i)=maxpna(maxat,i)+1
779 minpna(minat,i)=minpna(minat,i)+1
780 endif
781 endif
782
783 enddo
784
785 do j=3, taillr
786 nring(j,i) = tring(j)
787 enddo
788
789 !do j=3, TAILLR
790 ! write (6, '("s= ",i4,", j= ",i2,", nr(",i2,",",i4,")= ",i3)') i,j,j,i, NRING(j,i)
791 ! if (NRING(j,i) .gt. 0) then
792 ! do k=1, NRING(j,i)
793 ! write (6, *) " k= ",k,", R(o)= ",ORDRING(j,k,1:j)
794 ! enddo
795 ! endif
796 !enddo
797 j = rings_to_ogl(i, ari, nring, savring, ordring)
798
799#ifdef OPENMP
800 !$OMP ATOMIC
801#endif
802 ri = ri + j
803 002 continue
804
805enddo
806#ifdef OPENMP
807!$OMP END DO NOWAIT
808#endif
809
810001 continue
811
812if (alc) then
813 call show_error ("Impossible to allocate memory"//char(0), &
814 "Subroutine: KING_RING_SEARCH_STEPS"//char(0), "Table: "//alc_tab(1:len_trim(alc_tab))//char(0))
815endif
816
817if (allocated(tring)) deallocate (tring)
818if (allocated(the_ring)) deallocate (the_ring)
819if (allocated(savring)) deallocate (savring)
820if (allocated(ordring)) deallocate (ordring)
821if (allocated(cpat)) deallocate (cpat)
822if (allocated(vpat)) deallocate (vpat)
823if (allocated(rpat)) deallocate (rpat)
824if (allocated(res_list)) deallocate (res_list)
825if (allocated(indte)) deallocate (indte)
826if (allocated(indth)) deallocate (indth)
827if (allocated(apna)) deallocate (apna)
828
829#ifdef OPENMP
830!$OMP END PARALLEL
831#endif
832
833if (ri .eq. ns) ri = rings_to_ogl_menu(ari, nring)
834
835END SUBROUTINE
836
837INTEGER FUNCTION check_ring (THE_RING, FND, S_CR, RID, TAE, TAH, RSAVED, OSAVED, TRING, INDE, INDH, RESL)
838
839USE parameters
840
841IMPLICIT NONE
842
843TYPE (ring), DIMENSION(TAILLD), INTENT(IN) :: the_ring
844LOGICAL, INTENT(INOUT) :: fnd
845INTEGER, INTENT(IN) :: s_cr, rid
846INTEGER, INTENT(INOUT) :: tae, tah
847INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: rsaved, osaved
848INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: tring
849INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: inde, indh
850INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: resl
851INTEGER :: ix, iy, chaine
852LOGICAL :: is_ring
853LOGICAL :: homop
854DOUBLE PRECISION :: duv
855DOUBLE PRECISION, DIMENSION(3) :: rab, vab
856
857INTERFACE
858 DOUBLE PRECISION FUNCTION calcdij(R12, AT1, AT2, STEP_1, STEP_2, SID)
859 DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: r12
860 INTEGER, INTENT(IN) :: at1, at2, step_1, step_2, sid
861 END FUNCTION
862 SUBROUTINE save_this_ring (THE_RING, TLES, RSAVED, OSAVED, TRING, INDT, RESL)
863 USE parameters
864 TYPE (ring), DIMENSION(TAILLR), INTENT(IN) :: the_ring
865 INTEGER, INTENT(IN) :: tles
866 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: rsaved, osaved
867 INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: indt
868 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: resl, tring
869 END SUBROUTINE
870 SUBROUTINE del_this_ring (TLED, RSAVED, OSAVED, TRING, RESL, INDT)
871 USE parameters
872 INTEGER, INTENT(IN) :: tled
873 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: rsaved, osaved
874 INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: indt
875 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: resl, tring
876 END SUBROUTINE
877END INTERFACE
878
879 is_ring=.true.
880 homop=.false.
881 if (the_ring(rid)%ATOM.eq.the_ring(1)%ATOM .and. rid.ge.4) then
882
883 is_ring=.true.
884 homop=.false.
885 if (nsp .gt. 1) then
886 do ix=1, rid-1
887 if (the_ring(ix)%SPEC .eq. the_ring(ix+1)%SPEC) then
888 homop=.true.
889 exit
890 endif
891 enddo
892 endif
893
894 if (pbc) then
895
896 vab(:)=0.0d0
897 do ix=1, rid-1
898 if (ncells .gt. 1) then
899 duv = calcdij(rab,the_ring(ix)%ATOM,the_ring(ix+1)%ATOM,s_cr,s_cr,s_cr)
900 else
901 duv = calcdij(rab,the_ring(ix)%ATOM,the_ring(ix+1)%ATOM,s_cr,s_cr,1)
902 endif
903 do iy=1,3
904 vab(iy)=vab(iy)+rab(iy)
905 enddo
906 enddo
907 if (abs(vab(1)) .ge. 0.01 .or. abs(vab(2)) .ge. 0.01 .or. abs(vab(3)) .ge. 0.01) then
908 is_ring=.false.
909 endif
910
911 endif
912
913 if (is_ring) then
914
915 chaine = rid-1
916 if (allrings) then
917 if (.not.no_homo .or. .not.homop) then
918 fnd=.true.
919 call save_this_ring (the_ring, chaine, rsaved, osaved, tring, inde, resl)
920 endif
921 else
922 if (no_homo) then
923 if (chaine .le. tah) then
924 if (chaine.lt.tah .or. .not.homop) call del_this_ring (tah, rsaved, osaved, tring, resl, indh)
925 if (chaine.lt.tae) call del_this_ring (tae, rsaved, osaved, tring, resl, inde)
926 if (.not.homop .and. chaine.le.tae) then
927 fnd=.true.
928 tae=chaine
929 call save_this_ring (the_ring, tae, rsaved, osaved, tring, inde, resl)
930 else if (homop) then
931 if (chaine .lt. tae) then
932 fnd=.false.
933 tah=chaine
934 call save_this_ring (the_ring, tah, rsaved, osaved, tring, indh, resl)
935 endif
936 endif
937 endif
938 else
939 if (chaine .le. tae) then
940 if (chaine.lt.tae) call del_this_ring (tae, rsaved, osaved, tring, resl, inde)
941 fnd=.true.
942 tae=chaine
943 call save_this_ring (the_ring, tae, rsaved, osaved, tring, inde, resl)
944 endif
945 endif
946 endif
947 check_ring = 1
948
949 else
950
951 check_ring = 2
952
953 endif
954
955 else
956
957 check_ring = 0
958
959 endif
960
961END FUNCTION
962
963RECURSIVE SUBROUTINE inside_ring (THE_RING, FND, S_IR, AI_IR, RID, TAE, TAH, LRA, LRB, &
964 NRPAT, RSAVED, OSAVED, TRING, INDE, INDH, RESL, CPT, VPT)
965USE parameters
966
967IMPLICIT NONE
968
969TYPE (ring), DIMENSION(TAILLD), INTENT(INOUT) :: the_ring
970LOGICAL, INTENT(INOUT) :: fnd
971INTEGER, INTENT(IN) :: s_ir, ai_ir, rid
972INTEGER, INTENT(INOUT) :: tae, tah
973INTEGER, INTENT(IN) :: lra, lrb
974INTEGER, DIMENSION(NA), INTENT(INOUT) :: nrpat
975INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: rsaved, osaved
976INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: tring
977INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: inde, indh
978INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: resl
979INTEGER, DIMENSION(NA), INTENT(IN):: cpt
980INTEGER, DIMENSION(NA,MAXN), INTENT(IN) :: vpt
981INTEGER :: ind, res
982LOGICAL :: addsp
983
984INTERFACE
985 INTEGER FUNCTION check_ring (THE_RING, FND, S_CR, RID, TAE, TAH, RSAVED, OSAVED, TRING, INDE, INDH, RESL)
986 USE parameters
987 TYPE (ring), DIMENSION(TAILLR), INTENT(IN) :: the_ring
988 LOGICAL, INTENT(INOUT) :: fnd
989 INTEGER, INTENT(IN) :: s_cr, rid
990 INTEGER, INTENT(INOUT) :: tae, tah
991 INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: rsaved, osaved
992 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: tring
993 INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: inde, indh
994 INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: resl
995 END FUNCTION
996END INTERFACE
997
998if (rid-1 .lt. tae) then
999
1000 do while (the_ring(rid)%NEIGHBOR .ge. 1)
1001
1002 ind = vpt(ai_ir, the_ring(rid)%NEIGHBOR)
1003 if (nrpat(ind).eq.0) then
1004
1005 if (abab) then
1006 if (mod(rid,2).eq.0 .and. lot(ind).eq.lrb) then
1007 addsp=.true.
1008 elseif (mod(rid,2).ne.0 .and. lot(ind).eq.lra) then
1009 addsp=.true.
1010 else
1011 addsp=.false.
1012 endif
1013 else
1014 addsp=.true.
1015 endif
1016 if (addsp) then
1017
1018 the_ring(rid+1)%ATOM = ind
1019 the_ring(rid+1)%SPEC = lot(ind)
1020 the_ring(rid+1)%NEIGHBOR = cpt(ind)
1021 nrpat(ind)=1
1022 res = check_ring(the_ring, fnd, s_ir, rid+1, tae, tah, rsaved, osaved, tring, inde, indh, resl)
1023 if (tbr .or. alc) goto 001
1024 if (res .eq. 0) then
1025 call inside_ring (the_ring, fnd, s_ir, ind, rid+1, tae, tah, lra, lrb, &
1026 nrpat, rsaved, osaved, tring, inde, indh, resl, cpt, vpt)
1027 if (tbr .or. alc) goto 001
1028 endif
1029 nrpat(ind) = 0
1030
1031 endif
1032
1033 endif
1034
1035 the_ring(rid)%NEIGHBOR = the_ring(rid)%NEIGHBOR - 1
1036
1037 enddo
1038
1039endif
1040
1041001 continue
1042
1043END SUBROUTINE inside_ring
1044
1045SUBROUTINE save_this_ring (THE_RING, TLES, RSAVED, OSAVED, TRING, INDT, RESL)
1046
1047USE parameters
1048
1049IMPLICIT NONE
1050
1051TYPE (RING), DIMENSION(TAILLD), INTENT(IN) :: THE_RING
1052INTEGER, INTENT(IN) :: TLES
1053INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: RSAVED, OSAVED
1054INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: INDT
1055INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: RESL, TRING
1056INTEGER :: idx, idy
1057LOGICAL :: NEWRING
1058INTEGER, DIMENSION(TLES) :: TOTRI, TOSAV
1059
1060! A ring has been found, we need to check if it has already been found or not
1061
1062do idx=1, tles
1063
1064! Creation of two tab first for the numerical sorting of the list
1065! The other without sorting for output
1066 totri(idx)=the_ring(idx)%ATOM
1067 tosav(idx)=the_ring(idx)%ATOM
1068
1069enddo
1070
1071call tri(totri, tles)
1072
1073if (tring(tles) .ne. 0) then
1074
1075 do idx=1, tring(tles)
1076! do-loop on existing rings to check if the new on has already been found
1077
1078 newring=.false.
1079
1080 do idy=1, tles
1081 if (totri(idy) .ne. rsaved(tles,idx,idy)) then
1082 newring=.true.
1083 exit
1084 endif
1085 enddo
1086
1087 if (.not.newring) then
1088
1089 ! Already been found n-times, increment of the counter
1090 indt(idx)=indt(idx)+1
1091 exit
1092
1093 endif
1094
1095 enddo
1096
1097else
1098
1099 newring=.true.
1100
1101endif
1102
1103if (newring) then
1104
1105! A new ring has been found
1106! We save the data
1107 resl(tles)=resl(tles)+1
1108 tring(tles)=tring(tles)+1
1109 if (tring(tles) .gt. numa) then
1110 tbr=.true.
1111 goto 001
1112 endif
1113 indt(tring(tles))=indt(tring(tles))+1
1114 do idx=1, tles
1115 rsaved(tles,tring(tles),idx)=totri(idx)
1116 osaved(tles,tring(tles),idx)=tosav(idx)
1117 enddo
1118
1119endif
1120
1121001 continue
1122
1123END SUBROUTINE
1124
1125SUBROUTINE del_this_ring (TLED, RSAVED, OSAVED, TRING, RESL, INDT)
1126
1127USE parameters
1128
1129IMPLICIT NONE
1130
1131INTEGER, INTENT(IN) :: TLED
1132INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: RSAVED, OSAVED
1133INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: INDT
1134INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: RESL, TRING
1135INTEGER :: xdel, did
1136
1137! The max size of the ring possible for the triplet N1-At-N2 has down
1138! We have to delete the bigger rings already found for this triplet
1139if (resl(tled) .ge. 1) then
1140
1141 do xdel=0, resl(tled)-1
1142
1143 do did=1, taillr
1144
1145 rsaved(tled,tring(tled)-xdel,did)=0
1146 osaved(tled,tring(tled)-xdel,did)=0
1147
1148 enddo
1149
1150 do did=1, numa
1151
1152 indt(did)=0
1153
1154 enddo
1155
1156 enddo
1157
1158 tring(tled)=tring(tled)-resl(tled)
1159 resl(tled)=0
1160
1161endif
1162
1163END SUBROUTINE
min
#define min(a, b)
Definition global.h:81
max
#define max(a, b)
Definition global.h:80
show_error
void show_error(char *error, int val, GtkWidget *win)
show error message
Definition interface.c:293
parameters::allrings
logical allrings
Definition parameters.F90:111
parameters::contj
integer, dimension(:,:), allocatable contj
Definition parameters.F90:245
parameters::tbr
logical tbr
Definition parameters.F90:108
parameters::taillr
integer taillr
Definition parameters.F90:150
parameters::numa
integer numa
Definition parameters.F90:136
parameters::voisj
integer, dimension(:,:,:), allocatable voisj
Definition parameters.F90:294
parameters::vpat
integer, dimension(:,:), allocatable vpat
Definition parameters.F90:261
parameters::cpat
integer, dimension(:), allocatable cpat
Definition parameters.F90:215
recrings
integer function recrings(vid)
Definition resrings.F90:22
save_this_ring
subroutine save_this_ring(the_ring, tles, rsaved, osaved, tring, indt, resl)
Definition rings-king.F90:1046
king_rings
integer function king_rings()
Definition rings-king.F90:97
check_ring
integer function check_ring(the_ring, fnd, s_cr, rid, tae, tah, rsaved, osaved, tring, inde, indh, resl)
Definition rings-king.F90:838
del_this_ring
subroutine del_this_ring(tled, rsaved, osaved, tring, resl, indt)
Definition rings-king.F90:1126
king_ring_search_steps
subroutine king_ring_search_steps(ari)
Definition rings-king.F90:532
setup_cpat_vpat_ring
subroutine setup_cpat_vpat_ring(nat, str, cont, vois, cpt, vpt)
Definition rings-king.F90:22
inside_ring
recursive subroutine inside_ring(the_ring, fnd, s_ir, ai_ir, rid, tae, tah, lra, lrb, nrpat, rsaved, osaved, tring, inde, indh, resl, cpt, vpt)
Definition rings-king.F90:965
setup_cpat_vpat
subroutine setup_cpat_vpat(str, nat)
Definition rings-king.F90:62
rings_to_ogl
integer function rings_to_ogl(step, idsearch, nri, rsaved, osaved)
Definition rings_ogl.F90:22
rings_to_ogl_menu
integer function rings_to_ogl_menu(idsearch, nri)
Definition rings_ogl.F90:170
tri
subroutine tri(tab, dimtab)
Definition utils.F90:581
calcdij
double precision function calcdij(r12, at1, at2, step_1, step_2, sid)
Definition utils.F90:185