atomes 1.3.1
atomes: an atomic scale modeling tool box
Loading...
Searching...
No Matches
dmtx.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-2026 by CNRS and University of Strasbourg
15!
16!>
17!! @file dmtx.F90
18!! @short Distance matrix calculation
19!! @author Sébastien Le Roux <sebastien.leroux@ipcms.unistra.fr>
20
21SUBROUTINE print_this_pixel (TPS, TPIX, PIXL, pix, atoms)
22
23USE parameters
24
25IMPLICIT NONE
26
27TYPE (PIXEL), DIMENSION(TPS), INTENT(IN) :: TPIX
28TYPE (PIXEL), INTENT(IN) :: PIXL
29LOGICAL, INTENT(IN) :: atoms
30INTEGER, INTENT(IN) :: TPS, pix
31INTEGER :: PA, PB
32
33write (6, '("Pixel= ",i4)') pix
34do pa=1, pixl%NEIGHBOR
35 pb = pixl%IDNEIGH(pa)
36 if (atoms) then
37 write (6, '(" N(",i4,")= ",i10," contains: ",i4," atom(s)")') pa, pixl%IDNEIGH(pa), tpix(pb)%ATOMS
38 else
39 write (6, '(" N(",i4,")= ",i10)') pa, pixl%IDNEIGH(pa)
40 endif
41enddo
42write (6, *)
43if (atoms) then
44 do pa=1, pixl%ATOMS
45 write (6, '(" At(",i4,")= ",i10)') pa, pixl%ATOM_ID(pa)
46 enddo
47 write (6, *)
48endif
49
50END SUBROUTINE
51
52SUBROUTINE set_shift (shift, ai, bi, ci, npa, npb, npc, nab, abc)
53
54 IMPLICIT NONE
55
56 INTEGER, DIMENSION(3,3,3), INTENT(INOUT) :: shift
57 INTEGER, INTENT(IN) :: ai, bi, ci
58 INTEGER, INTENT(IN) :: npa, npb, npc
59 INTEGER, INTENT(IN) :: nab, abc
60
61 ! Standard pixel positions:
62
63 ! Lower level
64 !shift(1,1,1) = nab+npa+1
65 !shift(1,1,2) = nab+npa
66 !shift(1,1,3) = nab+npa-1
67 !shift(1,2,1) = nab+1
68 !shift(1,2,2) = nab
69 !shift(1,2,3) = nab-1
70 !shift(1,3,1) = nab-npa+1
71 !shift(1,3,2) = nab-npa
72 !shift(1,3,3) = nab-npa-1
73
74 ! Same level
75 !shift(2,1,1) = npa+1
76 !shift(2,1,2) = npa
77 !shift(2,1,3) = npa-1
78 !shift(2,2,1) = 1
79 !shift(2,2,2) = 0
80 !shift(2,2,3) = -1
81 !shift(2,3,1) = -npa+1
82 !shift(2,3,2) = -npa
83 !shift(2,3,3) = -npa-1
84
85 ! Higher level
86 !shift(3,1,1) = npa+1-nab
87 !shift(3,1,2) = npa-nab
88 !shift(3,1,3) = npa-1-nab
89 !shift(3,2,1) = 1-nab
90 !shift(3,2,2) = -nab
91 !shift(3,2,3) = -1-nab
92 !shift(3,3,1) = -npa+1-nab
93 !shift(3,3,2) = -npa-nab
94 !shift(3,3,3) = -npa-1-nab
95
96 ! The shift is to compensate PBC
97
98 if (ai .eq. 1) then
99 shift(1,:,:) = npa
100 else if (ai .eq. npa) then
101 shift(3,:,:) = - npa
102 endif
103
104 if (bi .eq. 1) then
105 shift(:,1,:) = shift(:,1,:) + nab
106 else if (bi .eq. npb) then
107 shift(:,3,:) = shift(:,3,:) - nab
108 endif
109
110 if (ci .eq. 1) then
111 shift(:,:,1) = shift(:,:,1) + abc
112 else if (ci .eq. npc) then
113 shift(:,:,3) = shift(:,:,3) - abc
114 endif
115
116END SUBROUTINE
117
118INTEGER FUNCTION getnbx(NP, NPS)
119
120USE parameters
121
122IMPLICIT NONE
123
124INTEGER, INTENT(IN) :: np, nps
125INTEGER :: pia, pib, pic, pid
126INTEGER :: tpixd, adapt_cut
127DOUBLE PRECISION :: icut, inid
128DOUBLE PRECISION :: mpsize
129DOUBLE PRECISION :: new_mpsize
130DOUBLE PRECISION :: targetdp
131DOUBLE PRECISION :: rhonum
132INTEGER :: memoid
133
134if (.not.pbc) then
135 do pia=1, 3
136 pmin(pia)=fullpos(1,pia,1)
137 pmax(pia)=pmin(pia)
138 enddo
139 do pia=1, nps
140 if (pia.eq.1) then
141 pid=2
142 else
143 pid=1
144 endif
145 do pib=pid, np
146 do pic=1, 3
147 pmin(pic) = min(fullpos(pib,pic,pia),pmin(pic))
148 pmax(pic) = max(fullpos(pib,pic,pia),pmax(pic))
149 enddo
150 enddo
151 enddo
152 do pia=1, 3
153 pmax(pia) = pmax(pia) - pmin(pia)
154 enddo
155else
156 pmax(:) = 0.0d0
157 do pib=1, ncells
158 do pia=1, 3
159 pmax(pia) = max(pmax(pia), the_box(pib)%modv(pia))
160 enddo
161 enddo
162 do pia=1, 3
163 pmin(pia) = - pmax(pia)/2.0
164 enddo
165endif
166
167#ifdef DEBUG
168 write (6, '("pmin:: x= ",f15.10,", y= ",f15.10,", z= ",f15.10)') pmin(1), pmin(2), pmin(3)
169 write (6, '("pmax:: x= ",f15.10,", y= ",f15.10,", z= ",f15.10)') pmax(1), pmax(2), pmax(3)
170#endif
171
172icut=0.0
173do pia=1, nsp
174 do pib=1, nsp
175 icut=max(icut,gr_tmp(pia,pib))
176 enddo
177enddo
178icut = sqrt(icut)
179icut = icut + 1.0d0
180! The pixel size must be larger than the cutoff
181! If the system is large with a low density the pixel box can be too big in number of pixels,
182! and it is not possible to allocate the required memory.
183! Therefore it is required to increase the cutoff slightly to reduce the number of pixels.
184! We need to look into the density of atoms per pixel, good target values range between 1.5 and 2.0.
185! With PBC this works providing that you do not have a large box that contains an isolated molecule,
186! That is an extremely low number density, thus we need to check that as well.
187targetdp=1.85d0
188mpsize=1.0d0
189! Always adjust cutoff if periodicity is not applied
190adapt_cut=2
191
192! Number density
193rhonum=np
194do pia=1, 3
195 rhonum=rhonum/pmax(pia)
196enddo
197if (rhonum.gt.0.01) adapt_cut=1
198
199do tpixd=1, adapt_cut
200
201! MPSIZE is the modifier
202 cutf=icut*mpsize
203 cutf=1.0d0/cutf
204 do pia=1, 3
205 isize(pia) = int(abs(pmax(pia))*cutf)
206 if (isize(pia) .lt. 3) isize(pia) = 1
207 enddo
208
209 getnbx = 1
210 if (calc_prings .and. pbc) then
211 getnbx = 5
212 do pia=1, 3
213 pmin(pia) = getnbx*pmin(pia)
214 pmax(pia) = getnbx*pmax(pia)
215 enddo
216 do pia=1, 3
217 isize(pia) = int(abs(pmax(pia))*cutf)
218 if (isize(pia) .lt. 3) isize(pia) = 1
219 enddo
220 endif
221
222 if (tpixd .eq. 1) then
223 pia = isize(1)*isize(2)*isize(3)
224 inid = np
225 inid = inid/pia
226 new_mpsize = (targetdp/inid)**(1.0/3.0)
227 if (new_mpsize .gt. mpsize) then
228 ! Adaptation only if the cutoff increases
229 ! Otherwise we would increase the number of pixels
230 mpsize = new_mpsize
231 endif
232 endif
233enddo
234
235#ifdef DEBUG
236
237 write (6, '("NBX= ",i10)') getnbx
238 write (6, '("isize:: x= ",I4,", y= ",I4,", z= ",I4)') isize(1), isize(2), isize(3)
239 memoid = guess_grid_memory_size(isize(1)*isize(2)*isize(3))
240 write (6, '("Estimated memory required to store the entire pixel grid in Mb : ",f15.10)') memoid/1000000.0d0
241 write (6, '("Estimated memory required to store the entire pixel grid in Mo : ",f15.10)') memoid/8000000.0d0
242
243 CONTAINS
244
245 INTEGER FUNCTION guess_grid_memory_size (abc)
246
247 INTEGER, INTENT(IN) :: abc
248
249 allocate(thepix(1))
250 allocate(thepix(1)%ATOM_ID(50))
251
252 guess_grid_memory_size = storage_size(thepix)*abc
253
254 deallocate(thepix(1)%ATOM_ID)
255 deallocate(thepix)
256
257 END FUNCTION
258
259#endif
260
261END FUNCTION
262
263#ifdef DEBUG
264SUBROUTINE print_pixel_grid ()
265
266USE parameters
267
268IMPLICIT NONE
269
270INTEGER :: pix, pid, cid, did, eid, ai, bi, ci
271INTEGER :: init_a, end_a
272INTEGER :: init_b, end_b
273INTEGER :: init_c, end_c
274INTEGER, DIMENSION(3) :: dim
275INTEGER, DIMENSION(3,3,3) :: shift
276LOGICAL :: boundary, keep_it
277
278dim(1) = -1
279dim(2) = 0
280dim(3) = 1
281call pix_info (isize(1), isize(2), isize(3))
282
283allocate(thepix(abc))
284
285pix=1
286
287do ci=1, isize(3)
288 do bi=1, isize(2)
289 do ai=1, isize(1)
290 shift(:,:,:) = 0
291 if (pbc .and. .not.calc_prings) then
292 call set_shift (shift, ai, bi, ci, isize(1), isize(2), isize(3), ab, abc)
293 boundary=.false.
294 else
295 if (ai.eq.1 .or. ai.eq.isize(1)) boundary=.true.
296 if (bi.eq.1 .or. bi.eq.isize(2)) boundary=.true.
297 if (ci.eq.1 .or. ci.eq.isize(3)) boundary=.true.
298 endif
299 init_a = 1
300 end_a = 3
301 if (isize(1) .eq. 1) then
302 init_a = 2
303 end_a = 2
304 endif
305 init_b = 1
306 end_b = 3
307 if (isize(2) .eq. 1) then
308 init_b = 2
309 end_b = 2
310 endif
311 init_c = 1
312 end_c = 3
313 if (isize(3) .eq. 1) then
314 init_c = 2
315 end_c = 2
316 endif
317 pid = 0
318 do cid=init_a, end_a
319 do did=init_b, end_b
320 do eid=init_c, end_c
321 keep_it=.true.
322 if (boundary) then
323 if (ai.eq.1 .and. cid.eq.1) then
324 keep_it=.false.
325 else if (ai.eq.isize(1) .and. cid.eq.3) then
326 keep_it=.false.
327 else if (bi.eq.1 .and. did.eq.1) then
328 keep_it=.false.
329 else if (bi.eq.isize(2) .and. did.eq.3) then
330 keep_it=.false.
331 else if (ci.eq.1 .and. eid.eq.1) then
332 keep_it=.false.
333 else if (ci.eq.isize(3) .and. eid.eq.3) then
334 keep_it=.false.
335 endif
336 endif
337 if (keep_it) then
338 pid = pid + 1
339 thepix(pix)%IDNEIGH(pid) = pix + dim(cid) + dim(did) * isize(1) + dim(eid) * ab
340 thepix(pix)%IDNEIGH(pid) = thepix(pix)%IDNEIGH(pid) + shift(cid,did,eid)
341 endif
342 enddo
343 enddo
344 enddo
345 thepix(pix)%NEIGHBOR = pid
346 call send_pix_info (pix-1, thepix(pix)%IDNEIGH, pid)
347 ! call PRINT_THIS_PIXEL (abc, THEPIX, THEPIX(pix), pix, .false.)
348 pix = pix + 1
349 enddo
350 enddo
351enddo
352
353END SUBROUTINE
354#endif
355
356LOGICAL FUNCTION distmtx(NAN, LAN, LOOKNGB, UPNGB)
357
358!
359! Compute the distance matrix - neighbors table
360!
361
362USE parameters
363USE mendeleiev
364
365#ifdef OPENMP
366!$ USE OMP_LIB
367#endif
368IMPLICIT NONE
369
370INTEGER, INTENT(IN) :: nan
371INTEGER, DIMENSION(NAN), INTENT(IN) :: lan
372LOGICAL, INTENT(IN) :: lookngb, upngb
373INTEGER :: sat
374INTEGER :: ra, rb, rc, rd, rf, rg, rh, ri, rj, rk, rl, rm
375INTEGER :: rn, ro, rp, rq, rs, rt, ru, rv, rw, rx, ry, rz
376INTEGER :: a_start, a_end
377INTEGER :: pix
378DOUBLE PRECISION :: dik
379DOUBLE PRECISION :: maxbd, minbd
380INTEGER, DIMENSION(:), ALLOCATABLE :: ba, bb
381INTEGER, DIMENSION(:), ALLOCATABLE :: ca, cb
382DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xc, yc, zc
383LOGICAL :: calcmat=.false.
384! Error message info !
385LOGICAL :: toom=.false.
386INTEGER :: tooi
387LOGICAL :: pixr=.false.
388INTEGER :: pout
389!
390LOGICAL :: is_clone, boundary, keep_it
391INTEGER, DIMENSION(3) :: pixpos
392DOUBLE PRECISION, DIMENSION(3) :: xyz, uvw
393INTEGER :: pid, cid, did, eid, fid, ai, bi, ci
394INTEGER :: init_a, end_a
395INTEGER :: init_b, end_b
396INTEGER :: init_c, end_c
397INTEGER, DIMENSION(3) :: dim = (/-1, 0, 1/)
398INTEGER, DIMENSION(3,3,3) :: shift
399#ifdef OPENMP
400INTEGER :: numth
401INTEGER :: thread_num
402INTEGER :: atom_start, atom_end
403INTEGER, DIMENSION(:), ALLOCATABLE :: atpix
404LOGICAL :: doatoms
405#endif
406
407INTERFACE
408 INTEGER FUNCTION getnbx (NP, NPS)
409 INTEGER, INTENT(IN) :: np, nps
410 END FUNCTION
411 LOGICAL FUNCTION dvtbox (ST, NAS, NAT, NPOS)
412 INTEGER, INTENT(IN) :: st, nas, nat
413 DOUBLE PRECISION, DIMENSION(NAT,3), INTENT(INOUT) :: npos
414 END FUNCTION
415#ifdef OPENMP
416 INTEGER FUNCTION get_thread_start (NOBJ, NTHREADS, THREAD_ID)
417 INTEGER, INTENT(IN) :: nobj, nthreads, thread_id
418 END FUNCTION
419 INTEGER FUNCTION get_thread_end (NOBJ, NTHREADS, THREAD_ID)
420 INTEGER, INTENT(IN) :: nobj, nthreads, thread_id
421 END FUNCTION
422#endif
423 DOUBLE PRECISION FUNCTION calcdij (R12, AT1, AT2, STEP_1, STEP_2, SID)
424 USE parameters
425 DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: r12
426 INTEGER, INTENT(IN) :: at1, at2, step_1, step_2, sid
427 END FUNCTION
428 DOUBLE PRECISION FUNCTION spheres_caps_volumes (DAB, RAP, RBP)
429 DOUBLE PRECISION, INTENT(IN) :: dab, rap, rbp
430 END FUNCTION
431END INTERFACE
432
433#ifdef DEBUG
434 write (6, '("In DMTX:: LOOKNGB= ",L1, ", NOHOM= ",L1," PBC= ",L1,", CUBIC= ",L1)') lookngb, nohp, pbc, overall_cubic
435#endif
436
437if (allocated(gr_tmp)) deallocate(gr_tmp)
438allocate(gr_tmp(nsp,nsp), stat=err)
439if (err .ne. 0) then
440 alc_tab="VOISJ"
441 alc=.true.
442 distmtx=.false.
443 goto 001
444endif
445
446do i=1, nsp
447 do j=1, nsp
448 if (gr_cut(i,j) .gt. gr_cutoff) then
449 gr_tmp(i,j)=gr_cutoff
450 else
451 gr_tmp(i,j)=gr_cut(i,j)
452 endif
453 if (.not.lookngb) gr_tmp(i,j) = gr_tmp(i,j) + 50.0
454 enddo
455enddo
456
457nbx = getnbx(nan, ns)
458
459if (pbc) then
460 nna = nbx**3*nan
461else
462 nna = nan
463endif
464
465if (lookngb) then
466 if (allocated(voisj)) deallocate(voisj)
467 allocate(voisj(maxn,nna,ns), stat=err)
468 if (err .ne. 0) then
469 alc_tab="VOISJ"
470 alc=.true.
471 distmtx=.false.
472 goto 001
473 endif
474 if (allocated(contj)) deallocate(contj)
475 allocate(contj(nna,ns), stat=err)
476 if (err .ne. 0) then
477 alc_tab="CONTJ"
478 alc=.true.
479 distmtx=.false.
480 goto 001
481 endif
482 contj(:,:)=0
483 voisj(:,:,:)=0
484endif
485
486if (lookngb) then
487 maxbd=0.0d0
488 minbd=10.0d0
489endif
490
491if (calc_prings) then
492 a_start = 1
493 a_end = nna
494else
495 a_start = nan *(nbx**3 - 1)/2 + 1
496 a_end = a_start + nan - 1
497endif
498
499ab = isize(1)*isize(2)
500abc = ab*isize(3)
501
502#ifdef OPENMP
503numth = omp_get_max_threads()
504doatoms=.false.
505if (ns.lt.numth) then
506 if (numth .ge. 2*(ns-1)) then
507 doatoms=.true.
508 else
509 numth=ns
510 endif
511endif
512
513if (all_atoms) doatoms=.true.
514
515toom=.false.
516
517#ifdef DEBUG
518 call print_pixel_grid ()
519#endif
520
521if (doatoms) then
522! OpemMP on Atoms
523 if (nan.lt.numth) numth=nan
524 distmtx=.true.
525 if (nbx.gt.1) then
526 allocate(poa(nna,3), stat=err)
527 if (err .ne. 0) then
528 alc_tab="POA"
529 alc=.true.
530 distmtx=.false.
531 goto 001
532 endif
533 endif
534
535 do sat=1, ns
536
537 if (nbx.gt.1) then
538 if (.not.dvtbox(sat, na, nna, poa)) then
539 distmtx=.false.
540 goto 001
541 endif
542 endif
543
544 if (allocated(thepix)) deallocate(thepix)
545 allocate(thepix(abc), stat=err)
546 if (err .ne. 0) then
547 alc_tab="THEPIX"
548 alc=.true.
549 distmtx = .false.
550 goto 001
551 endif
552
553 if (allocated(atpix)) deallocate(atpix)
554 allocate(atpix(nna), stat=err)
555 if (err .ne. 0) then
556 alc_tab="ATPIX"
557 alc=.true.
558 distmtx = .false.
559 goto 001
560 endif
561
562 ra = int(maxn*2.5)
563 ! 1) LOOKNGB: "RA" is the number of atom(s) in one pixel
564 ! with MAXN = 20, the '50' value seems high enough
565 ! to ensure to catch every atom in the pixel
566 ! this is ok providing that no funny cutoff is used
567 ! 2) .not.LOOKNGB: looking for tetrahedra, requires larger value:
568 if (.not.lookngb) ra = maxn*10
569 if (.not.lookngb .or. abc .le. 5) ra = min(ra*5, nna)
570 do rb=1, abc
571 thepix(rb)%ATOMS=0
572 thepix(rb)%TOCHECK=.false.
573 thepix(rb)%CHECKED=.false.
574 allocate(thepix(rb)%ATOM_ID(ra), stat=err)
575 if (err .ne. 0) then
576 alc_tab="THEPIX%ATOM_ID"
577 alc=.true.
578 distmtx = .false.
579 goto 001
580 endif
581 thepix(rb)%IDNEIGH(:) = 0
582 thepix(rb)%ATOM_ID(:) = 0
583 enddo
584
585 do ra=1, nna
586 if (.not. pbc) then
587 do rb=1, 3
588 if (isize(rb) .eq. 1) then
589 pixpos(rb) = 0
590 else
591 pixpos(rb) = int((fullpos(ra,rb,sat) - pmin(rb))*cutf)
592 if (pixpos(rb) .eq. isize(rb)) pixpos(rb) = pixpos(rb)-1
593 endif
594 enddo
595 else
596 if (nbx .gt. 1) then
597 if (ncells.gt.1) then
598 xyz = matmul(poa(ra,:), the_box(sat)%carttofrac/nbx)
599 else
600 xyz = matmul(poa(ra,:), the_box(1)%carttofrac/nbx)
601 endif
602 else
603 if (ncells.gt.1) then
604 xyz = matmul(fullpos(ra,:,sat), the_box(sat)%carttofrac)
605 else
606 xyz = matmul(fullpos(ra,:,sat), the_box(1)%carttofrac)
607 endif
608 endif
609 do rb=1, 3
610 uvw(rb) = xyz(rb) - floor(xyz(rb))
611 pixpos(rb) = int(uvw(rb)*isize(rb))
612 ! if (pixpos(RD) .eq. isize(RD)) pixpos(RD) = pixpos(RD)-1
613 enddo
614 endif
615 pix = pixpos(1) + pixpos(2)*isize(1) + pixpos(3)*ab + 1
616 if (pix.gt.abc .or. pix.le.0) then
617 write (6, '("Thread(id)= ",i4,", MD step= ",i8)') thread_num, sat
618 if (nbx .gt. 1) then
619 write (6, '("at= ",i7,", pos(x)= ",f15.10,", pos(y)= ",f15.10,", pos(z)= ",f15.10)') ra, poa(ra,:)
620 else
621 write (6, '("at= ",i7,", pos(x)= ",f15.10,", pos(y)= ",f15.10,", pos(z)= ",f15.10)') ra, fullpos(ra,:,sat)
622 endif
623 if (pbc) then
624 write (6, '("at= ",i7,", fra(x)= ",f15.10,", fra(y)= ",f15.10,", fra(z)= ",f15.10)') ra, xyz
625 write (6, '("at= ",i7,", cor(x)= ",f15.10,", cor(y)= ",f15.10,", cor(z)= ",f15.10)') ra, uvw
626 endif
627 write (6, '("at= ",i7,", pixpos(x)= ",i4,", pixpos(y)= ",i4,", pixpos(z)= ",i4)') ra, pixpos
628 write (6, '("at= ",i7,", pixpos= ",i10)') ra, pix
629 pixr=.true.
630 pout=pix
631 distmtx=.false.
632 goto 001
633 else
634 if ((pbc .and. ra.ge.a_start .and. ra.le.a_end) .or. .not.pbc) then
635 if (.not.thepix(pix)%TOCHECK) then
636 ai = pixpos(1) + 1
637 bi = pixpos(2) + 1
638 ci = pixpos(3) + 1
639 pid = 0
640 shift(:,:,:) = 0
641 boundary=.false.
642 if (pbc .and. .not.calc_prings) then
643 call set_shift (shift, ai, bi, ci, isize(1), isize(2), isize(3), ab, abc)
644 else
645 if (ai.eq.1 .or. ai.eq.isize(1)) boundary=.true.
646 if (bi.eq.1 .or. bi.eq.isize(2)) boundary=.true.
647 if (ci.eq.1 .or. ci.eq.isize(3)) boundary=.true.
648 endif
649 init_a = 1
650 end_a = 3
651 if (isize(1) .eq. 1) then
652 init_a = 2
653 end_a = 2
654 endif
655
656 init_b = 1
657 end_b = 3
658 if (isize(2) .eq. 1) then
659 init_b = 2
660 end_b = 2
661 endif
662
663 init_c = 1
664 end_c = 3
665 if (isize(3) .eq. 1) then
666 init_c = 2
667 end_c = 2
668 endif
669
670 do cid=init_a, end_a
671 do did=init_b, end_b
672 do eid=init_c, end_c
673 keep_it=.true.
674 if (boundary) then
675 if (ai.eq.1 .and. cid.eq.1) then
676 keep_it=.false.
677 else if (ai.eq.isize(1) .and. cid.eq.3) then
678 keep_it=.false.
679 else if (bi.eq.1 .and. did.eq.1) then
680 keep_it=.false.
681 else if (bi.eq.isize(2) .and. did.eq.3) then
682 keep_it=.false.
683 else if (ci.eq.1 .and. eid.eq.1) then
684 keep_it=.false.
685 else if (ci.eq.isize(3) .and. eid.eq.3) then
686 keep_it=.false.
687 endif
688 endif
689 if (keep_it) then
690 pid = pid+1
691 thepix(pix)%IDNEIGH(pid) = pix + dim(cid) + dim(did) * isize(1) + dim(eid) * ab + shift(cid,did,eid)
692 endif
693 enddo
694 enddo
695 enddo
696 thepix(pix)%NEIGHBOR = pid
697 endif
698 thepix(pix)%TOCHECK=.true.
699 endif
700 thepix(pix)%ATOMS = thepix(pix)%ATOMS + 1
701 thepix(pix)%ATOM_ID(thepix(pix)%ATOMS) = ra
702 endif
703 atpix(ra) = pix
704 enddo
705
706 ra = 0
707 rb = 0
708 !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
709 !$OMP& PRIVATE(THREAD_NUM, ATOM_START, ATOM_END, &
710 !$OMP& RC, RD, RF, RG, RH, RI, RJ, RK, RL, RM, &
711 !$OMP& RN, RO, RP, RQ, RS, RT, RU, RV, RW, RX, RY, RZ, ERR, ai, bi, ci, &
712 !$OMP& IS_CLONE, CALCMAT, Dij, Rij, Dik, BA, BB, CA, CB, XC, YC, ZC) &
713 !$OMP& SHARED(NUMTH, SAT, NS, NA, NNA, NAN, LAN, NSP, LOOKNGB, UPNGB, DISTMTX, &
714 !$OMP& NBX, PBC, THE_BOX, NCELLS, A_START, A_END, NOHP, MAXN, CUTF, &
715 !$OMP& POA, FULLPOS, Gr_TMP, CALC_PRINGS, MAXBD, MINBD, CONTJ, VOISJ, RA, RB, &
716 !$OMP& LA_COUNT, CORTA, CORNERA, EDGETA, EDGEA, DEFTA, DEFA, &
717 !$OMP& ALC, ALC_TAB, TOOM, TOOI, THEPIX, ATPIX)
718 thread_num = omp_get_thread_num()
719 atom_start = get_thread_start(nna, numth, thread_num)
720 atom_end = get_thread_end(nna, numth, thread_num)
721
722 do rc=atom_start, atom_end
723 rd = atpix(rc)
724 rf = rc - (rc/nan)*nan
725 if (rf .eq. 0) rf=nan
726 rg = lan(rf)
727 do rh=1, thepix(rd)%NEIGHBOR
728 ri = thepix(rd)%IDNEIGH(rh)
729 rj = thepix(ri)%ATOMS
730 do rk=1, rj
731 rl = thepix(ri)%ATOM_ID(rk)
732 if (rl .ne. rc) then
733 if ((rc.ge.a_start .and. rc.le.a_end) .or. (rl.ge.a_start .and. rl.le.a_end)) then
734 rm = rl - (rl/nan)*nan
735 if (rm .eq. 0) rm=nan
736 rn = lan(rm)
737 if (lookngb) then
738 if (.not.pbc .or. nbx.gt.1) then
739 if (((rc.ge.a_start .and. rc.le.a_end) .and. (rl.ge.a_start .and. rl.le.a_end)) .or. .not.pbc) then
740 is_clone=.false.
741 else
742 is_clone=.true.
743 endif
744 else
745 is_clone=.false.
746 endif
747 if (rg.ne.rn .or. .not.nohp) then
748 calcmat=.true.
749 else
750 calcmat=.false.
751 endif
752 if (calcmat) then
753 dij=0.0d0
754 if (nbx.gt.1) then
755 do rp=1,3
756 rij(rp) = poa(rc,rp) - poa(rl,rp)
757 dij=dij+rij(rp)**2
758 enddo
759 else
760 dik=0.0d0
761 do rp=1,3
762 rij(rp) = fullpos(rf,rp,sat) - fullpos(rm,rp,sat)
763 dik=dik+rij(rp)**2
764 enddo
765 if (ncells .gt. 1) then
766 dij = calcdij(rij,rf,rm,sat,sat,sat)
767 else
768 dij = calcdij(rij,rf,rm,sat,sat,1)
769 endif
770 if (dik-dij .gt.0.01d0) then
771 is_clone=.true.
772 endif
773 endif
774 if (dij .le. gr_tmp(rg,rn)) then
775 !$OMP ATOMIC
776 maxbd=max(dij,maxbd)
777 !$OMP ATOMIC
778 minbd=min(dij,minbd)
779 if (calc_prings) then
780 rp = rc
781 rq = rl
782 else
783 rp = rf
784 rq = rm
785 endif
786 contj(rp,sat)=contj(rp,sat)+1
787 if (contj(rp,sat) .gt. maxn) then
788 toom=.true.
789 tooi=rp
790 distmtx=.false.
791 goto 002
792 endif
793 voisj(contj(rp,sat),rp,sat)=rq
794 if (.not.calc_prings .and. upngb) then
795 if (is_clone) then
796 voisj(contj(rp,sat),rp,sat)=-rq
797 !$OMP ATOMIC
798 rb = rb + 1
799 else
800 !$OMP ATOMIC
801 ra = ra + 1
802 endif
803 endif
804 endif
805 endif
806 else
807 if (rg .eq. rn) then
808 do ro=1, nsp
809 if (la_count(rf,ro,sat).eq.4 .and. contj(rf,sat).eq.4) then
810 if (la_count(rm,ro,sat).eq.4 .and. contj(rm,sat).eq.4) then
811 rp=0
812 do rq=1, 4
813 rs = voisj(rq,rf,sat)
814 do rt=1, contj(rs,sat)
815 if (voisj(rt,rs,sat) .eq. rm) rp=rp+1
816 enddo
817 enddo
818 if (rp.eq.1) then
819 !$OMP ATOMIC
820 corta(rg,ro)=corta(rg,ro)+1
821 if (ns .gt. 1) then
822 !$OMP ATOMIC
823 cornera(rg,ro,sat)=cornera(rg,ro,sat)+1
824 endif
825 elseif (rp.eq.2) then
826 !$OMP ATOMIC
827 edgeta(rg,ro)=edgeta(rg,ro)+1
828 if (ns .gt. 1) then
829 !$OMP ATOMIC
830 edgea(rg,ro,sat)=edgea(rg,ro,sat)+1
831 endif
832 elseif (rp.ge.3) then
833 !$OMP ATOMIC
834 corta(rg,ro)=corta(rg,ro)+1
835 !$OMP ATOMIC
836 edgeta(rg,ro)=edgeta(rg,ro)+1
837 !$OMP ATOMIC
838 defta(rg,ro)=defta(rg,ro)+1
839 if (ns .gt. 1) then
840 !$OMP ATOMIC
841 cornera(rg,ro,sat)=cornera(rg,ro,sat)+1
842 !$OMP ATOMIC
843 edgea(rg,ro,sat)=edgea(rg,ro,sat)+1
844 !$OMP ATOMIC
845 defa(rg,ro,sat)=defa(rg,ro,sat)+1
846 endif
847 endif
848 endif
849 endif
850 enddo
851 endif
852 endif
853 endif
854 endif
855 enddo
856 enddo
857 enddo
858
859 002 continue
860 !$OMP END PARALLEL
861
862 if (allocated(thepix)) then
863 do rc=1, abc
864 if (allocated(thepix(rc)%ATOM_ID)) deallocate(thepix(rc)%ATOM_ID)
865 enddo
866 deallocate(thepix)
867 endif
868 if (allocated(atpix)) deallocate(atpix)
869 if (.not.distmtx) then
870 goto 001
871 endif
872 if (.not.lookngb) then
873 do rc=1, nsp
874 do rd=1, nsp
875 corta(rc,rd)=corta(rc,rd)/2
876 edgeta(rc,rd)=edgeta(rc,rd)/2
877 defta(rc,rd)=defta(rc,rd)/2
878 if (ns .gt. 1) then
879 cornera(rc,rd,sat)=cornera(rc,rd,sat)/2
880 edgea(rc,rd,sat)=edgea(rc,rd,sat)/2
881 defa(rc,rd,sat)=defa(rc,rd,sat)/2
882 endif
883 enddo
884 enddo
885 endif
886 if (.not.calc_prings .and.upngb) then
887 if (ra .gt. 0) then
888 if (allocated(ba)) deallocate(ba)
889 allocate(ba(ra), stat=err)
890 if (err .ne. 0) then
891 alc_tab="BA"
892 alc=.true.
893 distmtx=.false.
894 goto 001
895 endif
896 if (allocated(bb)) deallocate(bb)
897 allocate(bb(ra), stat=err)
898 if (err .ne. 0) then
899 alc_tab="BB"
900 alc=.true.
901 distmtx=.false.
902 goto 001
903 endif
904 endif
905 if (rb .gt. 0) then
906 if (allocated(ca)) deallocate(ca)
907 allocate(ca(rb), stat=err)
908 if (err .ne. 0) then
909 alc_tab="CA"
910 alc=.true.
911 distmtx=.false.
912 goto 001
913 endif
914 if (allocated(cb)) deallocate(cb)
915 allocate(cb(rb), stat=err)
916 if (err .ne. 0) then
917 alc_tab="CB"
918 alc=.true.
919 distmtx=.false.
920 goto 001
921 endif
922 if (allocated(xc)) deallocate(xc)
923 allocate(xc(rb), stat=err)
924 if (err .ne. 0) then
925 alc_tab="XC"
926 alc=.true.
927 distmtx=.false.
928 goto 001
929 endif
930 if (allocated(yc)) deallocate(yc)
931 allocate(yc(rb), stat=err)
932 if (err .ne. 0) then
933 alc_tab="YC"
934 alc=.true.
935 distmtx=.false.
936 goto 001
937 endif
938 if (allocated(zc)) deallocate(zc)
939 allocate(zc(rb), stat=err)
940 if (err .ne. 0) then
941 alc_tab="ZC"
942 alc=.true.
943 distmtx=.false.
944 goto 001
945 endif
946 endif
947 if (ra.gt.0 .or. rb.gt.0) then
948 ra=0
949 rb=0
950 do rc=1, nna
951 do rd=1, contj(rc,sat)
952 rf = voisj(rd,rc,sat)
953 if (rf .gt. 0) then
954 if (rf.gt.rc) then
955 ra=ra+1
956 ba(ra) = rc
957 bb(ra) = rf
958 endif
959 else
960 voisj(rd,rc,sat) = -rf
961 rf=-rf
962 if (rf.gt.rc) then
963 rb=rb+1
964 ca(rb) = rc
965 cb(rb) = rf
966 if (ncells .gt. 1) then
967 call calcrij (rc,rf,sat,sat,sat)
968 else
969 call calcrij (rc,rf,sat,sat,1)
970 endif
971 xc(rb) = rij(1)
972 yc(rb) = rij(2)
973 zc(rb) = rij(3)
974 endif
975 endif
976 enddo
977 enddo
978 endif
979 call update_bonds (0, sat-1, ra, ba, bb, xc, yc, zc)
980 call update_bonds (1, sat-1, rb, ca, cb, xc, yc, zc)
981 do rc=1, nna
982 call update_atom_neighbors (sat-1, rc-1, contj(rc,sat))
983 do rd=1, contj(rc,sat)
984 call update_this_neighbor (sat-1, rc-1, rd-1, voisj(rd,rc,sat))
985 enddo
986 enddo
987 endif
988
989 enddo ! En MD steps loop
990
991else
992
993 ! OpemMP on MD steps
994 distmtx=.true.
995 !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
996 !$OMP& PRIVATE(THEPIX, SAT, pix, pixpos, XYZ, UVW, shift, ai, bi, ci, &
997 !$OMP& RA, RB, RC, RD, RF, RG, RH, RI, RJ, RK, RL, RM, &
998 !$OMP& RN, RO, RP, RQ, RS, RT, RU, RV, RW, RX, RY, RZ, ERR, &
999 !$OMP& IS_CLONE, CALCMAT, Dij, Rij, Dik, BA, BB, CA, CB, XC, YC, ZC, POA, &
1000 !$OMP& pid, cid, did, eid, fid, init_a, end_a, init_b, end_b, init_c, end_c, boundary, keep_it) &
1001 !$OMP& SHARED(NUMTH, NS, NA, NNA, NAN, LAN, NSP, LOOKNGB, UPNGB, DISTMTX, &
1002 !$OMP& NBX, PBC, THE_BOX, NCELLS, isize, pmin, pmax, abc, ab, A_START, A_END, NOHP, MAXN, CUTF, &
1003 !$OMP& FULLPOS, CONTJ, VOISJ, Gr_TMP, CALC_PRINGS, MAXBD, MINBD, &
1004 !$OMP& LA_COUNT, CORTA, CORNERA, EDGETA, EDGEA, DEFTA, DEFA, &
1005 !$OMP& ALC, ALC_TAB, TOOM, TOOI, PIXR, POUT, dim)
1006#endif
1007
1008 if (allocated(thepix)) deallocate(thepix)
1009 allocate(thepix(abc), stat=err)
1010 if (err .ne. 0) then
1011 alc_tab="THEPIX"
1012 alc=.true.
1013 distmtx = .false.
1014#ifdef OPENMP
1015 goto 006
1016#else
1017 goto 001
1018#endif
1019 endif
1020
1021 ra = int(maxn*2.5)
1022 ! 1) LOOKNGB: "RA" is the number of atom(s) in one pixel
1023 ! with MAXN = 20, the '50' value seems high enough
1024 ! to ensure to catch every atom in the pixel
1025 ! this is ok providing that no funny cutoff is used
1026 ! 2) .not.LOOKNGB: looking for tetrahedra, requires larger value:
1027 if (.not.lookngb) ra = maxn*10
1028 if (.not.lookngb .or. abc .le. 5) ra = min(ra*5, nna)
1029 do rb=1, abc
1030 thepix(rb)%ATOMS=0
1031 thepix(rb)%TOCHECK=.false.
1032 thepix(rb)%CHECKED=.false.
1033 allocate(thepix(rb)%ATOM_ID(ra), stat=err)
1034 if (err .ne. 0) then
1035 alc_tab="THEPIX%ATOM_ID"
1036 alc=.true.
1037 distmtx = .false.
1038#ifdef OPENMP
1039 goto 006
1040#else
1041 goto 001
1042#endif
1043 endif
1044 thepix(rb)%IDNEIGH(:) = 0
1045 thepix(rb)%ATOM_ID(:) = 0
1046 enddo
1047#ifdef OPENMP
1048 !$OMP DO SCHEDULE(STATIC,NS/NUMTH)
1049 do sat=1, ns
1050
1051 if (.not.distmtx) goto 007
1052#else
1053 if (allocated(poa)) deallocate (poa)
1054 allocate(poa(nna,3), stat=err)
1055 if (err .ne. 0) then
1056 alc_tab="POA"
1057 alc=.true.
1058 distmtx=.false.
1059 goto 001
1060 endif
1061 do sat=1, ns
1062#endif
1063 if (nbx.gt.1) then
1064#ifdef OPENMP
1065 if (allocated(poa)) deallocate (poa)
1066 allocate(poa(nna,3), stat=err)
1067 if (err .ne. 0) then
1068 alc_tab="POA"
1069 alc=.true.
1070 distmtx=.false.
1071 goto 007
1072 endif
1073#endif
1074 if (.not.dvtbox(sat, na, nna, poa)) then
1075 distmtx=.false.
1076#ifdef OPENMP
1077 goto 007
1078#else
1079 goto 001
1080#endif
1081 endif
1082 endif
1083
1084 ! Clean THEPIX for each and every MD step
1085 do ra=1, abc
1086 thepix(ra)%ATOMS=0
1087 thepix(ra)%TOCHECK=.false.
1088 thepix(ra)%CHECKED=.false.
1089 thepix(ra)%IDNEIGH(:) = 0
1090 thepix(ra)%ATOM_ID(:) = 0
1091 enddo
1092 do ra=1, nna
1093 if (.not. pbc) then
1094 do rb=1, 3
1095 if (isize(rb) .eq. 1) then
1096 pixpos(rb) = 0
1097 else
1098 pixpos(rb) = int((fullpos(ra,rb,sat) - pmin(rb))*cutf)
1099 if (pixpos(rb) .eq. isize(rb)) pixpos(rb) = pixpos(rb)-1
1100 endif
1101 enddo
1102 else
1103 if (nbx .gt. 1) then
1104 if (ncells.gt.1) then
1105 xyz = matmul(poa(ra,:), the_box(sat)%carttofrac/nbx)
1106 else
1107 xyz = matmul(poa(ra,:), the_box(1)%carttofrac/nbx)
1108 endif
1109 else
1110 if (ncells.gt.1) then
1111 xyz = matmul(fullpos(ra,:,sat), the_box(sat)%carttofrac)
1112 else
1113 xyz = matmul(fullpos(ra,:,sat), the_box(1)%carttofrac)
1114 endif
1115 endif
1116 do rb=1, 3
1117 uvw(rb) = xyz(rb) - floor(xyz(rb))
1118 pixpos(rb) = int(uvw(rb)*isize(rb))
1119 ! if (pixpos(RB) .eq. isize(RB)) pixpos(RB) = pixpos(RB)-1
1120 enddo
1121 endif
1122 pix = pixpos(1) + pixpos(2)*isize(1) + pixpos(3)*ab + 1
1123 if (pix.gt.abc .or. pix.le.0) then
1124#ifdef OPENMP
1125 write (6, '("Thread(id)= ",i4,", MD step= ",i8)') omp_get_thread_num(), sat
1126#else
1127 write (6, '("MD step= ",i8)') sat
1128#endif
1129 if (nbx.gt.1) then
1130 write (6, '("at= ",i7,", pos(x)= ",f15.10,", pos(y)= ",f15.10,", pos(z)= ",f15.10)') ra, poa(ra,:)
1131 else
1132 write (6, '("at= ",i7,", pos(x)= ",f15.10,", pos(y)= ",f15.10,", pos(z)= ",f15.10)') ra, fullpos(ra,:,sat)
1133 endif
1134 if (pbc) then
1135 write (6, '("at= ",i7,", fra(x)= ",f15.10,", fra(y)= ",f15.10,", fra(z)= ",f15.10)') ra, xyz
1136 write (6, '("at= ",i7,", cor(x)= ",f15.10,", cor(y)= ",f15.10,", cor(z)= ",f15.10)') ra, uvw
1137 endif
1138 write (6, '("at= ",i7,", pixpos(x)= ",i4,", pixpos(y)= ",i4,", pixpos(z)= ",i4)') ra, pixpos
1139 write (6, '("at= ",i7,", pixpos= ",i10)') ra, pix
1140 pixr=.true.
1141 pout=pix
1142 distmtx=.false.
1143#ifdef OPENMP
1144 goto 007
1145#else
1146 goto 001
1147#endif
1148 else
1149 if ((pbc .and. ra.ge.a_start .and. ra.le.a_end) .or. .not.pbc) then
1150 if (.not.thepix(pix)%TOCHECK) then
1151 dim(1) = -1
1152 dim(2) = 0
1153 dim(3) = 1
1154 ai = pixpos(1) + 1
1155 bi = pixpos(2) + 1
1156 ci = pixpos(3) + 1
1157 pid = 0
1158 shift(:,:,:) = 0
1159 boundary=.false.
1160 if (pbc .and. .not.calc_prings) then
1161 call set_shift (shift, ai, bi, ci, isize(1), isize(2), isize(3), ab, abc)
1162 else
1163 if (ai.eq.1 .or. ai.eq.isize(1)) boundary=.true.
1164 if (bi.eq.1 .or. bi.eq.isize(2)) boundary=.true.
1165 if (ci.eq.1 .or. ci.eq.isize(3)) boundary=.true.
1166 endif
1167 init_a = 1
1168 end_a = 3
1169 if (isize(1) .eq. 1) then
1170 init_a = 2
1171 end_a = 2
1172 endif
1173
1174 init_b = 1
1175 end_b = 3
1176 if (isize(2) .eq. 1) then
1177 init_b = 2
1178 end_b = 2
1179 endif
1180
1181 init_c = 1
1182 end_c = 3
1183 if (isize(3) .eq. 1) then
1184 init_c = 2
1185 end_c = 2
1186 endif
1187
1188 do cid=init_a, end_a
1189 do did=init_b, end_b
1190 do eid=init_c, end_c
1191 keep_it=.true.
1192 if (boundary) then
1193 if (ai.eq.1 .and. cid.eq.1) then
1194 keep_it=.false.
1195 else if (ai.eq.isize(1) .and. cid.eq.3) then
1196 keep_it=.false.
1197 else if (bi.eq.1 .and. did.eq.1) then
1198 keep_it=.false.
1199 else if (bi.eq.isize(2) .and. did.eq.3) then
1200 keep_it=.false.
1201 else if (ci.eq.1 .and. eid.eq.1) then
1202 keep_it=.false.
1203 else if (ci.eq.isize(3) .and. eid.eq.3) then
1204 keep_it=.false.
1205 endif
1206 endif
1207 if (keep_it) then
1208 pid = pid+1
1209 thepix(pix)%IDNEIGH(pid) = pix + dim(cid) + dim(did) * isize(1) + dim(eid) * ab + shift(cid,did,eid)
1210 endif
1211 enddo
1212 enddo
1213 enddo
1214 thepix(pix)%NEIGHBOR = pid
1215 endif
1216 thepix(pix)%TOCHECK=.true.
1217 endif
1218 thepix(pix)%ATOMS = thepix(pix)%ATOMS + 1
1219 thepix(pix)%ATOM_ID(thepix(pix)%ATOMS) = ra
1220 endif
1221 enddo
1222
1223 ra = 0
1224 rb = 0
1225 do rc=1, abc
1226 if (thepix(rc)%TOCHECK) then
1227 rd = thepix(rc)%ATOMS
1228 if (rd .gt. 0) then
1229 do rf=1, thepix(rc)%NEIGHBOR
1230 rg = thepix(rc)%IDNEIGH(rf)
1231 if (.not.thepix(rg)%CHECKED) then
1232 rh = thepix(rg)%ATOMS
1233 if (rh .gt. 0) then
1234 if (rc .eq. rg) then
1235 ri=1
1236 else
1237 ri=0
1238 endif
1239 do rj=1, rd-ri
1240 rk = thepix(rc)%ATOM_ID(rj)
1241 do rm=ri*rj+1, rh
1242 rn = thepix(rg)%ATOM_ID(rm)
1243 if ((rk.ge.a_start .and. rk.le.a_end) .or. (rn.ge.a_start .and. rn.le.a_end)) then
1244 rp=rk - (rk/nan)*nan
1245 rq=rn - (rn/nan)*nan
1246 if (rp .eq. 0) rp=nan
1247 if (rq .eq. 0) rq=nan
1248 if (rp .ne. rq) then
1249 rl=lan(rp)
1250 ro=lan(rq)
1251 if (lookngb) then
1252 if (.not.pbc .or. nbx.gt.1) then
1253 if (((rk.ge.a_start .and. rk.le.a_end) .and. (rn.ge.a_start .and. rn.le.a_end)) .or. .not.pbc) then
1254 is_clone=.false.
1255 else
1256 is_clone=.true.
1257 endif
1258 else
1259 is_clone=.false.
1260 endif
1261 if (rl.ne.ro .or. .not.nohp) then
1262 calcmat=.true.
1263 else
1264 calcmat=.false.
1265 endif
1266 if (calcmat) then
1267 dij=0.0d0
1268 if (nbx.gt.1) then
1269 do rs=1,3
1270 rij(rs) = poa(rk,rs) - poa(rn,rs)
1271 dij=dij+rij(rs)**2
1272 enddo
1273 else
1274 dik=0.0d0
1275 do rs=1,3
1276 rij(rs) = fullpos(rp,rs,sat) - fullpos(rq,rs,sat)
1277 dik=dik+rij(rs)**2
1278 enddo
1279 if (ncells .gt. 1) then
1280 dij = calcdij(rij,rp,rq,sat,sat,sat)
1281 else
1282 dij = calcdij(rij,rp,rq,sat,sat,1)
1283 endif
1284 if (dik-dij .gt.0.01d0) then
1285 is_clone=.true.
1286 endif
1287 endif
1288 if (dij .le. gr_tmp(rl,ro)) then
1289#ifdef OPENMP
1290 !$OMP ATOMIC
1291#endif
1292 maxbd=max(dij,maxbd)
1293#ifdef OPENMP
1294 !$OMP ATOMIC
1295#endif
1296 minbd=min(dij,minbd)
1297 if (calc_prings) then
1298 rt = rk
1299 ru = rn
1300 else
1301 rt = rp
1302 ru = rq
1303 endif
1304
1305 contj(rt,sat)=contj(rt,sat)+1
1306 if (contj(rt,sat) .gt. maxn) then
1307 toom=.true.
1308 tooi=rt
1309 distmtx=.false.
1310#ifdef OPENMP
1311 goto 007
1312#else
1313 goto 001
1314#endif
1315 endif
1316 voisj(contj(rt,sat),rt,sat)=ru
1317 contj(ru,sat)=contj(ru,sat)+1
1318 if (contj(ru,sat) .gt. maxn) then
1319 toom=.true.
1320 tooi=ru
1321 distmtx=.false.
1322#ifdef OPENMP
1323 goto 007
1324#else
1325 goto 001
1326#endif
1327 endif
1328 voisj(contj(ru,sat),ru,sat)=rt
1329 if (.not.calc_prings .and. upngb) then
1330 if (is_clone) then
1331 rb = rb + 1
1332 voisj(contj(rt,sat),rt,sat)=-ru
1333 voisj(contj(ru,sat),ru,sat)=-rt
1334 else
1335 ra = ra + 1
1336 endif
1337 endif
1338 endif
1339 endif
1340 else
1341 if (rl .eq. ro) then
1342 do rv=1, nsp
1343 if (la_count(rp,rv,sat).eq.4 .and. contj(rp,sat).eq.4) then
1344 if (la_count(rq,rv,sat).eq.4 .and. contj(rq,sat).eq.4) then
1345 rw=0
1346 do rx=1, 4
1347 ry = voisj(rx,rp,sat)
1348 do rz=1, contj(ry,sat)
1349 if (voisj(rz,ry,sat) .eq. rq) rw=rw+1
1350 enddo
1351 enddo
1352 if (rw.eq.1) then
1353#ifdef OPENMP
1354 !$OMP ATOMIC
1355#endif
1356 corta(rl,rv)=corta(rl,rv)+1
1357 if (ns .gt. 1) cornera(rl,rv,sat)=cornera(rl,rv,sat)+1
1358 elseif (rw.eq.2) then
1359#ifdef OPENMP
1360 !$OMP ATOMIC
1361#endif
1362 edgeta(rl,rv)=edgeta(rl,rv)+1
1363 if (ns .gt. 1) edgea(rl,rv,sat)=edgea(rl,rv,sat)+1
1364 elseif (rw.ge.3) then
1365#ifdef OPENMP
1366 !$OMP ATOMIC
1367#endif
1368 corta(rl,rv)=corta(rl,rv)+1
1369#ifdef OPENMP
1370 !$OMP ATOMIC
1371#endif
1372 edgeta(rl,rv)=edgeta(rl,rv)+1
1373#ifdef OPENMP
1374 !$OMP ATOMIC
1375#endif
1376 defta(rl,rv)=defta(rl,rv)+1
1377 if (ns .gt. 1) then
1378 cornera(rl,rv,sat)=cornera(rl,rv,sat)+1
1379 edgea(rl,rv,sat)=edgea(rl,rv,sat)+1
1380 defa(rl,rv,sat)=defa(rl,rv,sat)+1
1381 endif
1382 endif
1383 endif
1384 endif
1385 enddo
1386 endif
1387 endif
1388 endif
1389 endif
1390 enddo
1391 enddo
1392 endif
1393 endif
1394 enddo
1395 thepix(rc)%CHECKED=.true.
1396 endif
1397 endif
1398 enddo
1399 if (.not.calc_prings .and.upngb) then
1400 if (ra .gt. 0) then
1401 if (allocated(ba)) deallocate(ba)
1402 allocate(ba(ra), stat=err)
1403 if (err .ne. 0) then
1404 alc_tab="BA"
1405 alc=.true.
1406 distmtx=.false.
1407#ifdef OPENMP
1408 goto 007
1409#else
1410 goto 001
1411#endif
1412 endif
1413 if (allocated(bb)) deallocate(bb)
1414 allocate(bb(ra), stat=err)
1415 if (err .ne. 0) then
1416 alc_tab="BB"
1417 alc=.true.
1418 distmtx=.false.
1419#ifdef OPENMP
1420 goto 007
1421#else
1422 goto 001
1423#endif
1424 endif
1425 endif
1426 if (rb .gt. 0) then
1427 if (allocated(ca)) deallocate(ca)
1428 allocate(ca(rb), stat=err)
1429 if (err .ne. 0) then
1430 alc_tab="CA"
1431 alc=.true.
1432 distmtx=.false.
1433#ifdef OPENMP
1434 goto 007
1435#else
1436 goto 001
1437#endif
1438 endif
1439 if (allocated(cb)) deallocate(cb)
1440 allocate(cb(rb), stat=err)
1441 if (err .ne. 0) then
1442 alc_tab="CB"
1443 alc=.true.
1444 distmtx=.false.
1445#ifdef OPENMP
1446 goto 007
1447#else
1448 goto 001
1449#endif
1450 endif
1451 if (allocated(xc)) deallocate(xc)
1452 allocate(xc(rb), stat=err)
1453 if (err .ne. 0) then
1454 alc_tab="XC"
1455 alc=.true.
1456 distmtx=.false.
1457#ifdef OPENMP
1458 goto 007
1459#else
1460 goto 001
1461#endif
1462 endif
1463 if (allocated(yc)) deallocate(yc)
1464 allocate(yc(rb), stat=err)
1465 if (err .ne. 0) then
1466 alc_tab="YC"
1467 alc=.true.
1468 distmtx=.false.
1469#ifdef OPENMP
1470 goto 007
1471#else
1472 goto 001
1473#endif
1474 endif
1475 if (allocated(zc)) deallocate(zc)
1476 allocate(zc(rb), stat=err)
1477 if (err .ne. 0) then
1478 alc_tab="ZC"
1479 alc=.true.
1480 distmtx=.false.
1481#ifdef OPENMP
1482 goto 007
1483#else
1484 goto 001
1485#endif
1486 endif
1487 endif
1488 if (ra.gt.0 .or. rb.gt.0) then
1489 ra=0
1490 rb=0
1491 do rc=1, nna
1492 do rd=1, contj(rc,sat)
1493 rf = voisj(rd,rc,sat)
1494 if (rf .gt. 0) then
1495 if (rf.gt.rc) then
1496 ra=ra+1
1497 ba(ra) = rc
1498 bb(ra) = rf
1499 endif
1500 else
1501 voisj(rd,rc,sat) = -rf
1502 rf=-rf
1503 if (rf.gt.rc) then
1504 rb=rb+1
1505 ca(rb) = rc
1506 cb(rb) = rf
1507 if (ncells .gt. 1) then
1508 dij = calcdij(rij,rc,rf,sat,sat,sat)
1509 else
1510 dij = calcdij(rij,rc,rf,sat,sat,1)
1511 endif
1512 xc(rb) = rij(1)
1513 yc(rb) = rij(2)
1514 zc(rb) = rij(3)
1515 endif
1516 endif
1517 enddo
1518 enddo
1519 endif
1520 call update_bonds (0, sat-1, ra, ba, bb, xc, yc, zc)
1521 call update_bonds (1, sat-1, rb, ca, cb, xc, yc, zc)
1522 do rc=1, nna
1523 call update_atom_neighbors (sat-1, rc-1, contj(rc,sat))
1524 do rd=1, contj(rc,sat)
1525 call update_this_neighbor (sat-1, rc-1, rd-1, voisj(rd,rc,sat))
1526 enddo
1527 enddo
1528 endif
1529
1530#ifdef OPENMP
1531 007 continue
1532#endif
1533
1534 enddo ! En MD steps loop
1535#ifdef OPENMP
1536 !$OMP END DO NOWAIT
1537
1538 006 continue
1539
1540 if (allocated(thepix)) then
1541 do rc=1, abc
1542 if (allocated(thepix(rc)%ATOM_ID)) deallocate(thepix(rc)%ATOM_ID)
1543 enddo
1544 deallocate(thepix)
1545 endif
1546 !$OMP END PARALLEL
1547
1548 if (.not.distmtx) goto 001
1549
1550endif
1551
1552#endif
1553
1554if (upngb) then
1555 maxbd = sqrt(maxbd)
1556 minbd = sqrt(minbd)
1557 if (maxbd-minbd .lt. 0.1) then
1558 minbd = minbd - 0.5;
1559 maxbd = maxbd + 0.5;
1560 endif
1561 call recup_dmin_dmax (minbd, maxbd)
1562endif
1563distmtx=.true.
1564
1565001 continue
1566
1567if (toom) call toomuch(tooi)
1568if (alc) then
1569 call show_error ("Impossible to allocate memory !"//char(0), &
1570 "Function: DMTX"//char(0), char(9)//"Table: "//alc_tab(1:len_trim(alc_tab))//char(0))
1571endif
1572if (pixr) call pixout (pout)
1573
1574if (allocated(thepix)) then
1575 do rc=1, abc
1576 if (allocated(thepix(rc)%ATOM_ID)) deallocate(thepix(rc)%ATOM_ID)
1577 enddo
1578 deallocate(thepix)
1579endif
1580if (allocated(poa)) deallocate(poa)
1581if (allocated(ba)) deallocate(ba)
1582if (allocated(bb)) deallocate(bb)
1583if (allocated(ca)) deallocate(ca)
1584if (allocated(cb)) deallocate(cb)
1585if (allocated(xc)) deallocate(xc)
1586if (allocated(yc)) deallocate(yc)
1587if (allocated(zc)) deallocate(zc)
1588
1589CONTAINS
1590
1591END FUNCTION
1592
1593SUBROUTINE pixout (PIX)
1594
1595 IMPLICIT NONE
1596
1597 INTEGER, INTENT(IN) :: PIX
1598 CHARACTER (LEN=15) :: IDPIX
1599 call charint(idpix, pix)
1600 call show_error ("Pixel size error !"//char(0), &
1601 char(9)//"Pixel out of bound(s): "//idpix(2:len_trim(idpix))//char(0), &
1602 "Function: DMTX"//char(0))
1603
1604END SUBROUTINE
1605
1606SUBROUTINE toomuch(ATI)
1607
1608IMPLICIT NONE
1609
1610INTEGER, INTENT(IN) :: ATI
1611CHARACTER (LEN=15) :: IDATI
1612
1613call charint(idati, ati)
1614
1615call show_error ("Too much neighbors for atom "//idati(2:len_trim(idati))//" (>20)"//char(0), &
1616 "Please check:"//achar(10)//char(9)//" - The lattice parameters"//achar(10) &
1617 //char(9)//" - The bond cutoff(s)"//achar(10)//char(9)//" - The atomic coordinates"//char(0), &
1618 "Function: DMTX"//char(0))
1619
1620END SUBROUTINE
1621
1622INTEGER (KIND=c_int) FUNCTION rundmtx (PRINGS, VNOHP, VUP) bind (C,NAME='rundmtx_')
1623
1624USE parameters
1625
1626IMPLICIT NONE
1627
1628INTEGER (KIND=c_int), INTENT(IN) :: prings, vnohp, vup
1629LOGICAL :: dmtxok=.false.
1630LOGICAL :: upng
1631
1632INTERFACE
1633 LOGICAL FUNCTION distmtx(NAN, LAN, LOOKNGB, UPNGB)
1634 INTEGER, INTENT(IN) :: nan
1635 INTEGER, DIMENSION(NAN), INTENT(IN) :: lan
1636 LOGICAL, INTENT(IN) :: lookngb, upngb
1637 END FUNCTION
1638END INTERFACE
1639
1640calc_prings=.false.
1641if (prings .eq. 1) calc_prings=.true.
1642nohp=.false.
1643if (vnohp .eq. 1) nohp=.true.
1644upng=.false.
1645if (vup .eq. 1) upng=.true.
1646
1647dmtxok = distmtx(na, lot, .true., upng)
1648calc_prings=.false.
1649
1650if (.not. dmtxok) then
1651 if (allocated(voisj)) deallocate(voisj)
1652 if (allocated(contj)) deallocate(contj)
1653 rundmtx=0
1654 goto 001
1655endif
1656
1657rundmtx=1
1658
1659001 continue
1660
1661END FUNCTION
1662
getnbx
integer function getnbx(np, nps)
Definition dmtx.F90:119
print_this_pixel
subroutine print_this_pixel(tps, tpix, pixl, pix, atoms)
Definition dmtx.F90:22
set_shift
subroutine set_shift(shift, ai, bi, ci, npa, npb, npc, nab, abc)
Definition dmtx.F90:53
rundmtx
integer(kind=c_int) function rundmtx(prings, vnohp, vup)
Definition dmtx.F90:1623
distmtx
logical function distmtx(nan, lan, lookngb, upngb)
Definition dmtx.F90:357
pixout
subroutine pixout(pix)
Definition dmtx.F90:1594
toomuch
subroutine toomuch(ati)
Definition dmtx.F90:1607
dvtbox
logical function dvtbox(st, nas, nat, npos)
Definition dvtb.F90:22
min
#define min(a, b)
Definition global.h:93
max
#define max(a, b)
Definition global.h:92
show_error
void show_error(char *error, int val, GtkWidget *win)
show error message
Definition interface.c:299
parameters::fullpos
double precision, dimension(:,:,:), allocatable fullpos
Definition parameters.F90:547
parameters::ncells
integer ncells
Definition parameters.F90:157
parameters::isize
integer, dimension(3) isize
Definition parameters.F90:196
parameters::gr_tmp
double precision, dimension(:,:), allocatable gr_tmp
Definition parameters.F90:515
parameters::contj
integer, dimension(:,:), allocatable contj
Definition parameters.F90:246
parameters::calc_prings
logical calc_prings
Definition parameters.F90:81
parameters::cutf
double precision cutf
Definition parameters.F90:413
parameters::pmax
double precision, dimension(3) pmax
Definition parameters.F90:447
parameters::thepix
type(pixel), dimension(:), allocatable thepix
Definition parameters.F90:653
parameters::the_box
type(lattice), dimension(:), allocatable, target the_box
Definition parameters.F90:668
parameters::nohp
logical nohp
Definition parameters.F90:69
parameters::pmin
double precision, dimension(3) pmin
Definition parameters.F90:447
parameters::voisj
integer, dimension(:,:,:), allocatable voisj
Definition parameters.F90:295
parameters::abc
integer abc
Definition parameters.F90:142
parameters::lot
integer, dimension(:), allocatable lot
Definition parameters.F90:187
parameters::pbc
logical pbc
Definition parameters.F90:67
parameters::nsp
integer nsp
Definition parameters.F90:147
get_thread_start
integer function get_thread_start(nobj, nthreads, thread_id)
Definition threads.F90:22
get_thread_end
integer function get_thread_end(nobj, nthreads, thread_id)
Definition threads.F90:37
charint
subroutine charint(word, num)
Definition utils.F90:32
calcdij
double precision function calcdij(r12, at1, at2, step_1, step_2, sid)
Definition utils.F90:185
calcrij
subroutine calcrij(at1, at2, step_1, step_2, sid)
Definition utils.F90:115