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