3 integer proba, conf, cm1, cm2, accept, nd
4 integer krok, i, n, j, g, ii1, ilosc, ijl
5 integer nnn, seed, proba1, proba2, proba3, proba0x,proba0y,proba0z
6 integer ilosc_atomow1, ilosc_atomow2, korekta1, korekta2
7 integer ile_reszt_mol1, ile_reszt_mol2, ir, peptide,numerr2(90000)
8 integer numera1(90000),numerr1(90000),numera2(90000)
9 integer glys1, glys2, glye1, glye2, ilee
10 real*8 rand, r, boxx, boxy, boxz, dist2acm, dist1acm, b(4), t(3)
11 real*8 minx1, miny1, minz1, maxx1, maxy1, maxz1
12 real*8 minx2, miny2, minz2, maxx2, maxy2, maxz2
13 real*8 distx1, disty1, distz1, distx2, disty2, distz2
14 real*8 xcm1, ycm1, zcm1, xcm2, ycm2, zcm2
15 real*8 dist1cm, dist1cmt, dist2cm, dist2cmt
16 real*8 cm2xr, cm2yr, cm2zr,x2r(90000),y2r(90000),z2r(90000),dist12
17 real*8 angle, c, s, x2new, y2new, z2new
18 character*60 plik1, plik2, plik3, plik4
20 real*8 x1(90000),y1(90000),z1(90000),x2(90000),y2(90000),z2(90000)
21 real*8 dd, wymiar, cmxx(900), cmyy(900), cmzz(900),xt2(90000,3)
24 character*3 r1s, r2s, r1e, r2e
25 character*5 atoma1(90000), atoma2(90000)
26 character*6 nazwaa1(90000), nazwar1(90000)
27 character*6 nazwaa2(90000), nazwar2(90000)
33 c write(*,*)"Sposob uruchamiania:"
34 c write(*,*)"nazwa_programu nazwa_pliku_inp_pdb1
35 c & nazwa_pliku_inp_pdb2 nazwa_pliku_out"
40 write(*,*)"program input_pdb1 input_pdb2
41 & [0-1] 0 for protein, 1 for peptide"
48 open(1,file=plik1,status="old")
49 open(2,file=plik2,status="old")
52 write(*,*) "All files read"
54 read(plik3(1:1),'(i1)') peptide
56 if (peptide.ne.0.and.peptide.ne.1) then
57 write(*,*) "Third parameter must be 0 or 1,
58 & where 0 is for a protein-protein docking,
59 & and 1 is for protein-peptide docking."
63 write(*,*) "Option", peptide, "selected
64 & (0 for protein, 1 for peptide)"
69 read(1,'(a54)',end=22) linia
70 if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then
71 write(4,'(a54)') linia
73 if (linia(1:3).eq.'TER') then
74 c write(2,'(a54)') linia
81 write(*,*) "Molecule 1 read"
86 read(2,'(a54)',end=23) linia
87 if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then
88 write(5,'(a54)') linia
90 if (linia(1:3).eq.'TER') then
91 c write(2,'(a54)') linia
97 write(*,*) "Molecule 2 read"
130 & atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),chain,numerr1(i),
134 c write(*,300) minx, maxx, miny, maxy, minz, maxz
136 if (x1(i).lt.minx1) minx1=x1(i)
137 if (y1(i).lt.miny1) miny1=y1(i)
138 if (z1(i).lt.minz1) minz1=z1(i)
140 if (x1(i).gt.maxx1) maxx1=x1(i)
141 if (y1(i).gt.maxy1) maxy1=y1(i)
142 if (z1(i).gt.maxz1) maxz1=z1(i)
144 c write(*,300) minx, maxx, miny, maxy, minz, maxz
146 c write(*,*) nazwaa1(i)
147 if (nazwaa1(i).eq.' CA ') then
151 ile_reszt_mol1 = ile_reszt_mol1 +1
154 if (atoma1(i).eq.'TER ') goto 111
155 if (atoma1(i).eq.'END ') goto 111
156 if (i.eq.1) write(3,'(a6)') "REMARK"
157 c write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),
158 c & numerr1(i),x1(i),y1(i),z1(i)
160 c write(*,*) i, x(i),y(i),z(i)
170 xcm1 = xcm1 / (ile_reszt_mol1*1.0)
171 ycm1 = ycm1 / (ile_reszt_mol1*1.0)
172 zcm1 = zcm1 / (ile_reszt_mol1*1.0)
185 if (nazwaa1(i).eq.' CA ') then
191 xcm1 = xcm1 / (ile_reszt_mol1*1.0)
192 ycm1 = ycm1 / (ile_reszt_mol1*1.0)
193 zcm1 = zcm1 / (ile_reszt_mol1*1.0)
198 c write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),
199 c & numerr1(i),x1(i)-xcm1,y1(i)-ycm1,z1(i)-zcm1
201 dist1cmt = dsqrt(((xcm1-x1(i))**2)+
205 c write(*,*) dist1cmt, dist1cm, i, cm1
206 if (dist1cm.gt.dist1cmt) then
215 dist1cmt = dsqrt(((x1(j)-x1(i))**2)+
216 & ((y1(j)-y1(i))**2)+
217 & ((z1(j)-z1(i))**2))
218 if (dist1acm.lt.dist1cmt) then
226 c write(3,304) "ENDMDL"
228 write(*,*) "Molecule 1 calculated"
239 & atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),chain,numerr2(i),
242 c write(*,300) minx, maxx, miny, maxy, minz, maxz
244 if (x2(i).lt.minx2) minx2=x2(i)
245 if (y2(i).lt.miny2) miny2=y2(i)
246 if (z2(i).lt.minz2) minz2=z2(i)
248 if (x2(i).gt.maxx2) maxx2=x2(i)
249 if (y2(i).gt.maxy2) maxy2=y2(i)
250 if (z2(i).gt.maxz2) maxz2=z2(i)
252 c write(*,300) minx, maxx, miny, maxy, minz, maxz
254 c write(*,*) nazwaa1(i)
255 if (nazwaa2(i).eq.' CA ') then
259 ile_reszt_mol2 = ile_reszt_mol2 + 1
262 if (atoma2(i).eq.'TER ') goto 111
263 if (atoma2(i).eq.'END ') goto 111
264 c if (i.eq.1) write(3,'(a6)') "REMARK"
265 c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),
266 c & numerr2(i),x2(i),y2(i),z2(i)
267 c write(*,*) i, x(i),y(i),z(i)
278 xcm2 = xcm2 / (ile_reszt_mol2*1.0)
279 ycm2 = ycm2 / (ile_reszt_mol2*1.0)
280 zcm2 = zcm2 / (ile_reszt_mol2*1.0)
286 c write(*,*) x2(i),y2(i),z2(i)
290 c write(*,*) xt2(i,1), xt2(i,2), xt2(i,3)
297 if (nazwaa2(i).eq.' CA ') then
303 xcm2 = xcm2 / (ile_reszt_mol2*1.0)
304 ycm2 = ycm2 / (ile_reszt_mol2*1.0)
305 zcm2 = zcm2 / (ile_reszt_mol2*1.0)
311 c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),
312 c & numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2
314 dist2cmt = dsqrt(((xcm2-x2(i))**2)+
318 c write(*,*) dist2cmt, dist2cm, i, cm2
319 if (dist2cm.gt.dist2cmt) then
329 dist2cmt = dsqrt(((x2(j)-x2(i))**2)+
330 & ((y2(j)-y2(i))**2)+
331 & ((z2(j)-z2(i))**2))
332 if (dist2acm.lt.dist2cmt) then
347 r1e=nazwar1(ile_reszt_mol1)
348 r2e=nazwar2(ile_reszt_mol2)
350 if (r1s.eq."GLY") then
352 write(*,*) "Chain 1 is starting from Gly ", r1s, glys1
355 if (r2s.eq."GLY") then
357 write(*,*) "Chain 2 is starting from Gly ", r2s, glys2
360 if (r1e.eq."GLY") then
362 write(*,*) "Chain 1 is ending at Gly ", r1e, glye1
366 if (r2e.eq."GLY") then
368 write(*,*) "Chain 2 is ending from Gly ", r2e, glye2
374 if(glys1.eq.1) ilee=ilee-1
375 if(glys2.eq.1) ilee=ilee-1
376 if(glye1.eq.1) ilee=ilee-1
377 if(glye2.eq.1) ilee=ilee-1
378 write(*,*) "How many dummy atoms:", ilee
382 write(*,*) "Molecule 2 calculated"
383 write(*,*) "Number of atoms:", ilosc_atomow1, ilosc_atomow2
385 write(*,*) "Min and max coordinates in the systems"
386 write(*,300) minx1, maxx1, miny1, maxy1, minz1, maxz1
387 write(*,300) minx2, maxx2, miny2, maxy2, minz2, maxz2
389 write(*,*) "Distances in axes x, y, z"
390 write(*,301) distx1, disty1, distz1
391 write(*,301) distx2, disty2, distz2
393 boxx = (distx1 + distx2)*1.5 + 20.0
394 boxy = (disty1 + disty2)*1.5 + 20.0
395 boxz = (distz1 + distz2)*1.5 + 20.0
397 c write(*,*) "Boxsize [x, y, z]"
398 c write(*,301) boxx, boxy, boxz
400 write(*,*) "Maximum distances:"
401 write(*,*) dist1acm, dist2acm
403 write(*,*) "Boxsize (suggested)"
404 write(*,*) (dist1acm+dist2acm)*1.2+20
405 open(61,file='boxsize.txt')
406 write(61,302) ((dist1acm+dist2acm)*1.2+20)
409 write(*,*) "Center of masses:"
410 write(*,301) xcm1, ycm1, zcm1
411 write(*,301) xcm2, ycm2, zcm2
413 write(*,*) "Center of masses atoms:"
414 write(*,800) atoma1(cm1),numera1(cm1),nazwaa1(cm1),nazwar1(cm1),
415 & " ",numerr1(cm1),x1(cm1),y1(cm1),z1(cm1)
416 write(*,800) atoma2(cm2),numera2(cm2),nazwaa2(cm2),nazwar2(cm2),
417 & " ",numerr2(cm2),x2(cm2),y2(cm2),z2(cm2)
420 write(*,*) "Number of amino-acid residues in two proteins:"
421 write(*,*) ile_reszt_mol1, ile_reszt_mol2
424 c ---------------------------------------------------------------------
426 c ---------------------------------------------------------------------
428 korekta1 = numerr1(1)
429 write(*,*) "Original number of the first residue in
430 & the first protein:"
431 write(*,*) numerr1(1)
433 write(*,*) "Correction of the first residue in the first protein:"
437 c wymiar=(distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2)
441 c write(*,*) "Minimum required distance between centers of mass:"
444 boxx = (7+distx1 + dist2acm)/2.0
445 boxy = (7+disty1 + dist2acm)/2.0
446 boxz = (7+distz1 + dist2acm)/2.0
447 write(*,*) "Rotation space (+/-):", boxx, boxy, boxz
463 write (plik3, "(A6,I1,a4)") "model0", ii1, ".pdb"
465 write (plik3, "(A5,I2,a4)") "model", ii1, ".pdb"
468 write(*,*) "Conformation will be saved to file:", plik3
470 c plik3a = plik3//ilosc
475 write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),"A",
476 & (numerr1(i)-korekta1+1),
480 c write(3,304) "ENDMDL"
482 if (peptide.eq.1) then
483 wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/4
486 if (peptide.eq.0) then
487 wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/7
491 write(*,*) "Initial minimum distance between centers of mass:",
495 c ROTATIONS OF THE SECOND MOLECULE
500 call ranorND(b,nd,ilosc_atomow2,ii1)
504 call rot_trans(xt2,t,b,xf,ilosc_atomow2)
514 if (conf.gt.1) goto 133
520 seed = (60*time_a(2)+time_a(3))*proba
524 c write(*,*)rand(0),rand(0),rand(0)
525 c write(*,*)rand(0),rand(0),rand(0)
527 cm2xr = (rand(0)*2*boxx)-(boxx)
528 cm2yr = (rand(0)*2*boxy)-(boxy)
529 cm2zr = (rand(0)*2*boxz)-(boxz)
530 c write(*,*) cm2xr, cm2yr, cm2zr
531 c write(*,300)rand(0),rand(0),rand(0)
533 c write(*,*) ABS(cm2xr+xcm1), (dist1acm+dist2acm+3)
535 c if (ABS(cm2xr+xcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
540 c if (ABS(cm2yr+ycm1).gt.((dist1acm+dist2acm+5)*0.8)) then
545 c if (ABS(cm2zr+zcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
558 c SPRAWDZANIE CZY NIE JEST ZBYT PODOBNE DO JUZ WYLOSOWANEGO
566 dd=sqrt((cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+
567 & (cm2zr-cmzz(ir))**2)
569 c write(*,*) "Center of mass is", dd, "A away from", ir
570 c write(*,*) "Minimum distance required:", wymiar, "A"
572 c write(*,*) ii1, ir, proba3, dd, wymiar
574 if (dd.lt.wymiar) then
576 if (proba3.gt.1.and.peptide.eq.1) wymiar=wymiar*0.99992
577 if (proba3.gt.1.and.peptide.eq.0) wymiar=wymiar*0.99992
578 c write(*,*) proba3, wymiar
588 c SPRAWDZANIE ODLEGLOSCI
591 c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),
592 c & numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2
595 dist12 = sqrt(((x2r(i)-x1(j))**2)+
596 & ((y2r(i)-y1(j))**2)+
597 & ((z2r(i)-z1(j))**2))
598 c write(*,*) dist12, cm2xr, cm2yr, cm2zr
599 if (dist12.lt.2.5) then
604 if (dist12.le.3.5) accept=1
605 c if (dist12.le.3.5) accept=1
606 c if (dist12.gt.3.5) proba2=proba2+1
610 c write(*,*) "No overlaps", proba
611 if (accept.eq.0) goto 200
612 c write(*,*) "No overlaps, contact", conf, proba
613 c write(3,305) "MODEL", conf
618 write(*,*) "Conformation", ii1, "Is accepted after",
619 & proba0x+proba0y+proba0z,
620 & proba1, proba2, proba3, proba, "attempts. ",
621 & "Dist between CoM of previous=", wymiar
624 write(*,*)"Extr",ABS(cm2xr+xcm1), ABS(cm2yr+ycm1), ABS(cm2zr+zcm1)
634 korekta2 = numerr2(1)-1
635 c write(*,*) "Original number of the first residue in
636 c & the second protein:"
637 c write(*,*) numerr2(1)
639 c write(*,*)"Correction of the first residue in the second protein:"
640 c write(*,*) korekta2
642 write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),"B",
643 & (numerr2(i)+ile_reszt_mol1-korekta2),
644 & x2r(i),y2r(i),z2r(i)
647 c write(3,304) "ENDMDL"
648 c write(3,305) "MODEL", conf
662 write (plik4, "(A9,I1,a4)") "template0", ii1, ".sco"
664 write (plik4, "(A8,I2,a4)") "template", ii1, ".sco"
670 if (peptide.eq.1) then
671 do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee)
673 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
674 if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
675 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
676 if (ijl.ge.(ile_reszt_mol1+2))
677 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
681 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
682 if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
683 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
684 if (ijl.ge.(ile_reszt_mol1+1))
685 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
691 if (peptide.eq.0) then
692 do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee)
695 if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.0) then
696 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
697 if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
698 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
699 if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3))
700 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
701 if (ijl.gt.(ile_reszt_mol1+3).and.
702 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+3))
703 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
704 if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+4))
705 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
710 if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.0) then
711 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
712 if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
713 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
714 if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+2))
715 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
716 if (ijl.gt.(ile_reszt_mol1+2).and.
717 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2))
718 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
719 if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3))
720 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
724 if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.0) then
725 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
726 if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
727 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
728 if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2))
729 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
730 if (ijl.gt.(ile_reszt_mol1+2).and.
731 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2))
732 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
733 if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3))
734 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
737 c chain 1 first Gly and chain 1 last Gly
738 if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.0) then
739 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
740 if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
741 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
742 if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1))
743 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
744 if (ijl.gt.(ile_reszt_mol1+1).and.
745 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1))
746 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
747 if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2))
748 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
752 if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.1) then
753 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
754 if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
755 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
756 if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2))
757 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
758 if (ijl.gt.(ile_reszt_mol1+2).and.
759 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2))
760 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
761 if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3))
762 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
765 c chain 1 first Gly chain 2 first Gly
766 if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.1) then
767 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
768 if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
769 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
770 if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1))
771 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
772 if (ijl.gt.(ile_reszt_mol1+1).and.
773 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1))
774 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
775 if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2))
776 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
779 c chain 1 last Gly chain 2 first Gly
780 if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.1) then
781 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
782 if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
783 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
784 c if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3))
785 c & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
786 if (ijl.gt.(ile_reszt_mol1+1).and.
787 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1))
788 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
789 if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2))
790 & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
794 if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.1) then
795 if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
796 if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
797 & write(7,*) ijl, "1.0 1.0 1.0 ", "1"
798 c if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3))
799 c & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
800 if (ijl.gt.(ile_reszt_mol1).and.
801 & ijl.le.(ile_reszt_mol2+ile_reszt_mol1))
802 & write(7,*) ijl, "1.0 1.0 1.0 ", "2"
803 c if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+1))
804 c & write(7,*) ijl, "1.0 1.0 1.0 ", "0"
815 c koniec petli po ilosci struktur
819 c write(*,*) i, cmxx(i), cmyy(i), cmzz(i)
820 c dd=(cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+(cm2zr-cmyy(ir))**2
825 c dd=sqrt((cmxx(i)-cmxx(j))**2+(cmyy(i)-cmyy(j))**2+
826 c & (cmzz(i)-cmzz(j))**2)
835 800 format(a4,2X,i5,1X,a4,1X,a3,1X,a1,i4,4X,3f8.3)
844 SUBROUTINE SRAND(ISEED)
846 C This subroutine sets the integer seed to be used with the
847 C companion RAND function to the value of ISEED. A flag is
848 C set to indicate that the sequence of pseudo-random numbers
849 C for the specified seed should start from the beginning.
851 COMMON /SEED/JSEED,IFRST
860 C This function returns a pseudo-random number for each invocation.
861 C It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal
862 C standard number generator whose Pascal code appears in the article:
864 C Park, Steven K. and Miller, Keith W., "Random Number Generators:
865 C Good Ones are Hard to Find", Communications of the ACM,
868 PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773,
871 COMMON /SEED/JSEED,IFRST
872 INTEGER HVLUE, LVLUE, TESTV, NEXTN
875 IF (IFRST .EQ. 0) THEN
880 HVLUE = NEXTN / MOBYMP
881 LVLUE = MOD(NEXTN, MOBYMP)
882 TESTV = MPLIER*LVLUE - MOMDMP*HVLUE
883 IF (TESTV .GT. 0) THEN
886 NEXTN = TESTV + MODLUS
888 RAND = REAL(NEXTN)/REAL(MODLUS)
893 COMMON /SEED/JSEED,IFRST
895 DATA JSEED,IFRST/123456789,0/
899 Subroutine ranorND(b,nd,ilosc_atomow2,ii1)
901 C generates a random vector on a unit sphere onC an ND-dimensional space.
902 C algorithm 40 page 410.
903 C the algorithm assumes that the formula for bz is valid, i.e
904 C bz=(1.0-2.0*ransq),
905 C and also that, ransq can be computed as as the sum of the independent random
906 C numbers to the second power (expression 3).
907 dimension ran(nd-1), b(nd)
909 c write(*,*) "In Subroutine ranorND(b,nd), nd=", nd
910 c write(*,*) "ilosc_atomow2=", ilosc_atomow2
911 c do j=1,ilosc_atomow2
915 do while (ransq.ge.1)
919 ran(i)=1.0-2.0*rand(0)
920 c write(*,*) i, ran(i)
921 rsq=ran(i)*ran(i)+rsq ! (3)
925 c write(*,*) j, ransq
927 ranh=2.0*sqrt(1.0-ransq)
928 c write(*,*) j, ranh, ransq
931 c write(*,*) ran(i), ranh
933 b(nd)=(1.0-2.0*ransq)
934 c write(*,*) b(1,j), b(2,j), b(3,j), b(4,j)
939 subroutine rot_trans(xt2,t,q,xf,ilosc_atomow2)
940 C x coordinates, t translation vector, q quaternion
941 double precision xt2(90000,3),xf(90000,3)
942 double precision t(3),q(4)
944 write(*,*) "Rotation vector:", q(1), q(2), q(3), q(4)
945 r11= q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4)
946 r22= q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4)
947 r33= q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4)
948 r12= 2*(q(2)*q(3)-q(1)*q(4))
949 r21= 2*(q(2)*q(3)+q(1)*q(4))
950 r13= 2*(q(2)*q(4)+q(1)*q(3))
951 r31= 2*(q(2)*q(4)-q(1)*q(3))
952 r23= 2*(q(3)*q(4)-q(1)*q(2))
953 r32= 2*(q(3)*q(4)+q(1)*q(2))
955 c write(*,*) r11, r22, r33, r12, r21, r13, r31, r23, r32
958 xf(i,1)= r11*xt2(i,1)+r12*xt2(i,2)+r13*xt2(i,3) + t(1)
959 xf(i,2)= r21*xt2(i,1)+r22*xt2(i,2)+r23*xt2(i,3) + t(2)
960 xf(i,3)= r31*xt2(i,1)+r32*xt2(i,2)+r33*xt2(i,3) + t(3)
961 c write(*,*) i, xt2(i,1), xt2(i,2), xt2(i,3)
962 c write(*,*) i, xf(i,1), xf(i,2), xf(i,3)