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