program generator implicit none integer proba, conf, cm1, cm2, accept, nd integer krok, i, n, j, g, ii1, ilosc, ijl integer nnn, seed, proba1, proba2, proba3, proba0x,proba0y,proba0z integer ilosc_atomow1, ilosc_atomow2, korekta1, korekta2 integer ile_reszt_mol1, ile_reszt_mol2, ir, peptide,numerr2(90000) integer numera1(90000),numerr1(90000),numera2(90000) integer glys1, glys2, glye1, glye2, ilee integer no_models, seed1, seed11, seed12, seed13 real*8 rand, r, boxx, boxy, boxz, dist2acm, dist1acm, b(4), t(3) real*8 minx1, miny1, minz1, maxx1, maxy1, maxz1 real*8 minx2, miny2, minz2, maxx2, maxy2, maxz2 real*8 distx1, disty1, distz1, distx2, disty2, distz2 real*8 xcm1, ycm1, zcm1, xcm2, ycm2, zcm2 real*8 dist1cm, dist1cmt, dist2cm, dist2cmt real*8 cm2xr, cm2yr, cm2zr,x2r(90000),y2r(90000),z2r(90000),dist12 real*8 angle, c, s, x2new, y2new, z2new character*60 plik1, plik2, plik3, plik4, arg4, arg5 character*54 linia real*8 x1(90000),y1(90000),z1(90000),x2(90000),y2(90000),z2(90000) real*8 dd, wymiar, cmxx(900), cmyy(900), cmzz(900),xt2(90000,3) real*8 xf(90000,3) character*1 chain character*3 r1s, r2s, r1e, r2e character*5 atoma1(90000), atoma2(90000) character*6 nazwaa1(90000), nazwar1(90000) character*6 nazwaa2(90000), nazwar2(90000) korekta1 = 1 korekta2 = 1 c write(*,*)"Sposob uruchamiania:" c write(*,*)"nazwa_programu nazwa_pliku_inp_pdb1 c & nazwa_pliku_inp_pdb2 nazwa_pliku_out" c write(*,*)"" nnn = iargc() if (nnn.lt.5) then write(*,*)"program input_pdb1 input_pdb2 & [0-1] number_of_models seed_number & 0 for protein, 1 for peptide" stop endif call getarg(1,plik1) call getarg(2,plik2) call getarg(3,plik3) call getarg(4,arg4) call getarg(5,arg5) open(1,file=plik1,status="old") open(2,file=plik2,status="old") open(4,file='temp1') open(5,file='temp2') write(*,*) "All files read" read(plik3(1:1),'(i1)') peptide read(arg4,'(i2)') no_models read(arg5,'(i19)') seed1 if (peptide.ne.0.and.peptide.ne.1) then write(*,*) "Third parameter must be 0 or 1, & where 0 is for a protein-protein docking, & and 1 is for protein-peptide docking." STOP endif write(*,*) "Option", peptide, "selected & (0 for protein, 1 for peptide)" c MOLEKULA 1 do while (.true.) read(1,'(a54)',end=22) linia if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then write(4,'(a54)') linia endif if (linia(1:3).eq.'TER') then c write(2,'(a54)') linia goto 22 endif enddo 22 continue close(4) write(*,*) "Molecule 1 read" c MOLEKULA 2 do while (.true.) read(2,'(a54)',end=23) linia if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then write(5,'(a54)') linia endif if (linia(1:3).eq.'TER') then c write(2,'(a54)') linia goto 23 endif enddo 23 continue close(5) write(*,*) "Molecule 2 read" open(4,file='temp1') open(5,file='temp2') minx1=999.99 miny1=999.99 minz1=999.99 maxx1=-999.99 maxy1=-999.99 maxz1=-999.99 minx2=999.99 miny2=999.99 minz2=999.99 maxx2=-999.99 maxy2=-999.99 maxz2=-999.99 c MOLEKULA 1 i=1 xcm1=0 ycm1=0 zcm1=0 ile_reszt_mol1=0 do while (.true.) read(4,800,end=111) & atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),chain,numerr1(i), & x1(i),y1(i),z1(i) c write(*,300) minx, maxx, miny, maxy, minz, maxz if (x1(i).lt.minx1) minx1=x1(i) if (y1(i).lt.miny1) miny1=y1(i) if (z1(i).lt.minz1) minz1=z1(i) if (x1(i).gt.maxx1) maxx1=x1(i) if (y1(i).gt.maxy1) maxy1=y1(i) if (z1(i).gt.maxz1) maxz1=z1(i) c write(*,300) minx, maxx, miny, maxy, minz, maxz c write(*,*) nazwaa1(i) if (nazwaa1(i).eq.' CA ') then xcm1 = xcm1 + x1(i) ycm1 = ycm1 + y1(i) zcm1 = zcm1 + z1(i) ile_reszt_mol1 = ile_reszt_mol1 +1 endif if (atoma1(i).eq.'TER ') goto 111 if (atoma1(i).eq.'END ') goto 111 if (i.eq.1) write(3,'(a6)') "REMARK" c write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i), c & numerr1(i),x1(i),y1(i),z1(i) c write(*,*) i, x(i),y(i),z(i) i=i+1 enddo 111 continue distx1=maxx1-minx1 disty1=maxy1-miny1 distz1=maxz1-minz1 ilosc_atomow1 = i-1 xcm1 = xcm1 / (ile_reszt_mol1*1.0) ycm1 = ycm1 / (ile_reszt_mol1*1.0) zcm1 = zcm1 / (ile_reszt_mol1*1.0) do i=1,ilosc_atomow1 x1(i)=x1(i)-xcm1 y1(i)=y1(i)-ycm1 z1(i)=z1(i)-zcm1 enddo xcm1=0 ycm1=0 zcm1=0 do i=1,ilosc_atomow1 if (nazwaa1(i).eq.' CA ') then xcm1 = xcm1 + x1(i) ycm1 = ycm1 + y1(i) zcm1 = zcm1 + z1(i) endif enddo xcm1 = xcm1 / (ile_reszt_mol1*1.0) ycm1 = ycm1 / (ile_reszt_mol1*1.0) zcm1 = zcm1 / (ile_reszt_mol1*1.0) dist1cm=999.9 dist1acm=0.0 do i=1,ilosc_atomow1 c write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i), c & numerr1(i),x1(i)-xcm1,y1(i)-ycm1,z1(i)-zcm1 dist1cmt = dsqrt(((xcm1-x1(i))**2)+ & ((ycm1-y1(i))**2)+ & ((zcm1-z1(i))**2)) c write(*,*) dist1cmt, dist1cm, i, cm1 if (dist1cm.gt.dist1cmt) then dist1cm=dist1cmt cm1=i endif enddo do i=1,ilosc_atomow1 do j=i,ilosc_atomow1 dist1cmt = dsqrt(((x1(j)-x1(i))**2)+ & ((y1(j)-y1(i))**2)+ & ((z1(j)-z1(i))**2)) if (dist1acm.lt.dist1cmt) then dist1acm=dist1cmt endif enddo enddo c write(3,303) "TER" c write(3,304) "ENDMDL" write(*,*) "Molecule 1 calculated" c MOLEKULA 2 i=1 xcm2=0 ycm2=0 zcm2=0 ile_reszt_mol2=0 do while (.true.) read(5,800,end=112) & atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),chain,numerr2(i), & x2(i),y2(i),z2(i) c write(*,300) minx, maxx, miny, maxy, minz, maxz if (x2(i).lt.minx2) minx2=x2(i) if (y2(i).lt.miny2) miny2=y2(i) if (z2(i).lt.minz2) minz2=z2(i) if (x2(i).gt.maxx2) maxx2=x2(i) if (y2(i).gt.maxy2) maxy2=y2(i) if (z2(i).gt.maxz2) maxz2=z2(i) c write(*,300) minx, maxx, miny, maxy, minz, maxz c write(*,*) nazwaa1(i) if (nazwaa2(i).eq.' CA ') then xcm2 = xcm2 + x2(i) ycm2 = ycm2 + y2(i) zcm2 = zcm2 + z2(i) ile_reszt_mol2 = ile_reszt_mol2 + 1 endif if (atoma2(i).eq.'TER ') goto 111 if (atoma2(i).eq.'END ') goto 111 c if (i.eq.1) write(3,'(a6)') "REMARK" c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i), c & numerr2(i),x2(i),y2(i),z2(i) c write(*,*) i, x(i),y(i),z(i) i=i+1 enddo 112 continue distx2=maxx2-minx2 disty2=maxy2-miny2 distz2=maxz2-minz2 ilosc_atomow2 = i-1 xcm2 = xcm2 / (ile_reszt_mol2*1.0) ycm2 = ycm2 / (ile_reszt_mol2*1.0) zcm2 = zcm2 / (ile_reszt_mol2*1.0) do i=1,ilosc_atomow2 x2(i)=x2(i)-xcm2 y2(i)=y2(i)-ycm2 z2(i)=z2(i)-zcm2 c write(*,*) x2(i),y2(i),z2(i) xt2(i,1)=x2(i) xt2(i,2)=y2(i) xt2(i,3)=z2(i) c write(*,*) xt2(i,1), xt2(i,2), xt2(i,3) enddo xcm2=0 ycm2=0 zcm2=0 do i=1,ilosc_atomow2 if (nazwaa2(i).eq.' CA ') then xcm2 = xcm2 + x2(i) ycm2 = ycm2 + y2(i) zcm2 = zcm2 + z2(i) endif enddo xcm2 = xcm2 / (ile_reszt_mol2*1.0) ycm2 = ycm2 / (ile_reszt_mol2*1.0) zcm2 = zcm2 / (ile_reszt_mol2*1.0) dist2cm=999.9 dist2acm=0.0 do i=1,ilosc_atomow2 c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i), c & numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2 dist2cmt = dsqrt(((xcm2-x2(i))**2)+ & ((ycm2-y2(i))**2)+ & ((zcm2-z2(i))**2)) c write(*,*) dist2cmt, dist2cm, i, cm2 if (dist2cm.gt.dist2cmt) then dist2cm=dist2cmt cm2=i endif enddo do i=1,ilosc_atomow2 do j=i,ilosc_atomow2 dist2cmt = dsqrt(((x2(j)-x2(i))**2)+ & ((y2(j)-y2(i))**2)+ & ((z2(j)-z2(i))**2)) if (dist2acm.lt.dist2cmt) then dist2acm=dist2cmt endif enddo enddo c write(3,303) "TER" glys1=0 glys2=0 glye1=0 glye2=0 r1s=nazwar1(1) r2s=nazwar2(1) r1e=nazwar1(ile_reszt_mol1) r2e=nazwar2(ile_reszt_mol2) if (r1s.eq."GLY") then glys1=1 write(*,*) "Chain 1 is starting from Gly ", r1s, glys1 endif if (r2s.eq."GLY") then glys2=1 write(*,*) "Chain 2 is starting from Gly ", r2s, glys2 endif if (r1e.eq."GLY") then glye1=1 write(*,*) "Chain 1 is ending at Gly ", r1e, glye1 endif if (r2e.eq."GLY") then glye2=1 write(*,*) "Chain 2 is ending from Gly ", r2e, glye2 endif c checking for Gly ilee=4 if(glys1.eq.1) ilee=ilee-1 if(glys2.eq.1) ilee=ilee-1 if(glye1.eq.1) ilee=ilee-1 if(glye2.eq.1) ilee=ilee-1 write(*,*) "How many dummy atoms:", ilee write(*,*) "" write(*,*) "Molecule 2 calculated" write(*,*) "Number of atoms:", ilosc_atomow1, ilosc_atomow2 write(*,*) "Min and max coordinates in the systems" write(*,300) minx1, maxx1, miny1, maxy1, minz1, maxz1 write(*,300) minx2, maxx2, miny2, maxy2, minz2, maxz2 write(*,*) "Distances in axes x, y, z" write(*,301) distx1, disty1, distz1 write(*,301) distx2, disty2, distz2 boxx = (distx1 + distx2)*1.5 + 20.0 boxy = (disty1 + disty2)*1.5 + 20.0 boxz = (distz1 + distz2)*1.5 + 20.0 c write(*,*) "Boxsize [x, y, z]" c write(*,301) boxx, boxy, boxz write(*,*) "Maximum distances:" write(*,*) dist1acm, dist2acm write(*,*) "Boxsize (suggested)" write(*,*) (dist1acm+dist2acm)*1.2+20 open(61,file='boxsize.txt') write(61,302) ((dist1acm+dist2acm)*1.2+20) close (61) write(*,*) "Center of masses:" write(*,301) xcm1, ycm1, zcm1 write(*,301) xcm2, ycm2, zcm2 write(*,*) "Center of masses atoms:" write(*,800) atoma1(cm1),numera1(cm1),nazwaa1(cm1),nazwar1(cm1), & " ",numerr1(cm1),x1(cm1),y1(cm1),z1(cm1) write(*,800) atoma2(cm2),numera2(cm2),nazwaa2(cm2),nazwar2(cm2), & " ",numerr2(cm2),x2(cm2),y2(cm2),z2(cm2) write(*,*) "Number of amino-acid residues in two proteins:" write(*,*) ile_reszt_mol1, ile_reszt_mol2 c --------------------------------------------------------------------- c KONIEC WCZYTYWANIA c --------------------------------------------------------------------- korekta1 = numerr1(1) write(*,*) "Original number of the first residue in & the first protein:" write(*,*) numerr1(1) write(*,*) "Correction of the first residue in the first protein:" write(*,*) korekta1 c wymiar=(distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2) c wymiar=wymiar/12.0 c wymiar=10.0 c write(*,*) "Minimum required distance between centers of mass:" c write(*,*) wymiar boxx = (6.5+distx1 + dist2acm)/2.0 boxy = (6.5+disty1 + dist2acm)/2.0 boxz = (6.5+distz1 + dist2acm)/2.0 write(*,*) "Rotation space (+/-):", boxx, boxy, boxz do ii1=1,no_models proba = 1 proba0x = 1 proba0y = 1 proba0z = 1 proba1 = 1 proba2 = 1 proba3 = 1 conf = 1 if (ii1.lt.10) then write (plik3, "(A6,I1,a4)") "model0", ii1, ".pdb" else write (plik3, "(A5,I2,a4)") "model", ii1, ".pdb" endif write(*,*) "Conformation will be saved to file:", plik3 c plik3a = plik3//ilosc open(3,file=plik3) do i=1,ilosc_atomow1 write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),"A", & (numerr1(i)-korekta1+1), & x1(i),y1(i),z1(i) enddo write(3,303) "TER" c write(3,304) "ENDMDL" if (peptide.eq.1) then wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/4 endif if (peptide.eq.0) then wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/8 endif if (ii1.eq.1) then write(*,*) "Initial minimum distance between centers of mass:", & wymiar endif c ROTATIONS OF THE SECOND MOLECULE nd=4 t(1)=0.0 t(2)=0.0 t(3)=0.0 call ranorND(b,nd,ilosc_atomow2,ii1,seed1) c do i=1,4 c write(*,*) b(i) c enddo call rot_trans(xt2,t,b,xf,ilosc_atomow2) do i=1,ilosc_atomow2 x2(i)=xf(i,1) y2(i)=xf(i,2) z2(i)=xf(i,3) enddo c END OF ROTATIONS do while (.true.) if (conf.gt.1) goto 133 seed = ABS((seed1+proba)/(ii1**2)) call srand(seed) r = rand(seed) seed11=ABS(seed+ii1-proba)*proba seed12=ABS(seed+ii1-(proba**2))*proba seed13=ABS(seed+ii1-(proba))*proba**2 c write(*,*) seed11, seed12, seed13 cm2xr = (rand(seed11)*2*boxx)-(boxx) cm2yr = (rand(seed12)*2*boxy)-(boxy) cm2zr = (rand(seed13)*2*boxz)-(boxz) do i=1,ilosc_atomow2 x2r(i)=x2(i)-cm2xr y2r(i)=y2(i)-cm2yr z2r(i)=z2(i)-cm2zr enddo c SPRAWDZANIE CZY NIE JEST ZBYT PODOBNE DO JUZ WYLOSOWANEGO cmxx(ii1)=cm2xr cmyy(ii1)=cm2yr cmzz(ii1)=cm2zr if (ii1.gt.1) then do ir=1,ii1-1 dd=sqrt((cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+ & (cm2zr-cmzz(ir))**2) c write(*,*) "Center of mass is", dd, "A away from", ir c write(*,*) "Minimum distance required:", wymiar, "A" c write(*,*) ii1, ir, proba3, dd, wymiar if (dd.lt.wymiar) then proba3=proba3+1 if (proba3.gt.1.and.peptide.eq.1) wymiar=wymiar*0.99992 if (proba3.gt.1.and.peptide.eq.0) wymiar=wymiar*0.99992 c write(*,*) proba3, wymiar goto 200 endif enddo endif c SPRAWDZANIE ODLEGLOSCI accept=0 do i=1,ilosc_atomow2 c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i), c & numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2 do j=1,ilosc_atomow1 dist12 = sqrt(((x2r(i)-x1(j))**2)+ & ((y2r(i)-y1(j))**2)+ & ((z2r(i)-z1(j))**2)) c write(*,*) dist12, cm2xr, cm2yr, cm2zr if (dist12.lt.2.5) then proba1=proba1+1 goto 200 endif if (dist12.le.3.5) accept=1 c if (dist12.le.3.5) accept=1 c if (dist12.gt.3.5) proba2=proba2+1 enddo enddo c write(*,*) "No overlaps", proba if (accept.eq.0) goto 200 c write(*,*) "No overlaps, contact", conf, proba c write(3,305) "MODEL", conf c conf = conf+1 write(*,*) "Conformation", ii1, "Is accepted after", & proba0x+proba0y+proba0z, & proba1, proba2, proba3, proba, "attempts. ", & "Dist between CoM of previous=", wymiar c Dodatkowe info write(*,*)"Extr",ABS(cm2xr+xcm1), ABS(cm2yr+ycm1), ABS(cm2zr+zcm1) conf=conf+1 c write (*,*) korekta2 = numerr2(1)-1 c write(*,*) "Original number of the first residue in c & the second protein:" c write(*,*) numerr2(1) c write(*,*)"Correction of the first residue in the second protein:" c write(*,*) korekta2 do i=1,ilosc_atomow2 write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),"B", & (numerr2(i)+ile_reszt_mol1-korekta2), & x2r(i),y2r(i),z2r(i) enddo write(3,303) "TER" c write(3,304) "ENDMDL" c write(3,305) "MODEL", conf write(3,303) "END" 200 continue proba = proba + 1 enddo 133 continue close (3) c call sleep(1) if (ii1.lt.10) then write (plik4, "(A9,I1,a4)") "template0", ii1, ".sco" else write (plik4, "(A8,I2,a4)") "template", ii1, ".sco" endif open (7,file=plik4) if (peptide.eq.1) then do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee) if (glys1.eq.0) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif if (glys1.eq.1) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+1)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif enddo endif if (peptide.eq.0) then do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee) c write(*,*) "ijl" c no Gly if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.0) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1+3).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+3)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+4)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif c chain 1 first Gly if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.0) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1+2).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif c chain 1 last Gly if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.0) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1+2).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif c chain 1 first Gly and chain 1 last Gly if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.0) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1+1).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif c chain 2 first Gly if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.1) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1+2).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif c chain 1 first Gly chain 2 first Gly if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.1) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1+1).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif c chain 1 last Gly chain 2 first Gly if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.1) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" c if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3)) c & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1+1).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2)) & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif c all Gly if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.1) then if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) & write(7,*) ijl, "1.0 1.0 1.0 ", "1" c if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3)) c & write(7,*) ijl, "1.0 1.0 1.0 ", "0" if (ijl.gt.(ile_reszt_mol1).and. & ijl.le.(ile_reszt_mol2+ile_reszt_mol1)) & write(7,*) ijl, "1.0 1.0 1.0 ", "2" c if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+1)) c & write(7,*) ijl, "1.0 1.0 1.0 ", "0" endif enddo endif close (7) c koniec petli po ilosci struktur enddo c do i=1,20 c write(*,*) i, cmxx(i), cmyy(i), cmzz(i) c dd=(cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+(cm2zr-cmyy(ir))**2 c enddo c do i=1,20 c do j=1,20 c dd=sqrt((cmxx(i)-cmxx(j))**2+(cmyy(i)-cmyy(j))**2+ c & (cmzz(i)-cmzz(j))**2) c write(*,*) i,j,dd c enddo c enddo 300 format(6f8.2) 301 format(3f8.2) 302 format(f8.2) 800 format(a4,2X,i5,1X,a4,1X,a3,1X,a1,i4,4X,3f8.3) 303 format(a3) 304 format(a6) 305 format(a5,1X,i8) end SUBROUTINE SRAND(ISEED) C C This subroutine sets the integer seed to be used with the C companion RAND function to the value of ISEED. A flag is C set to indicate that the sequence of pseudo-random numbers C for the specified seed should start from the beginning. C COMMON /SEED/JSEED,IFRST C JSEED = ISEED IFRST = 0 C RETURN END REAL FUNCTION RAND() C C This function returns a pseudo-random number for each invocation. C It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal C standard number generator whose Pascal code appears in the article: C C Park, Steven K. and Miller, Keith W., "Random Number Generators: C Good Ones are Hard to Find", Communications of the ACM, C October, 1988. C PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773, + MOMDMP=2836) C COMMON /SEED/JSEED,IFRST INTEGER HVLUE, LVLUE, TESTV, NEXTN SAVE NEXTN C IF (IFRST .EQ. 0) THEN NEXTN = JSEED IFRST = 1 ENDIF C HVLUE = NEXTN / MOBYMP LVLUE = MOD(NEXTN, MOBYMP) TESTV = MPLIER*LVLUE - MOMDMP*HVLUE IF (TESTV .GT. 0) THEN NEXTN = TESTV ELSE NEXTN = TESTV + MODLUS ENDIF RAND = REAL(NEXTN)/REAL(MODLUS) C RETURN END BLOCKDATA RANDBD COMMON /SEED/JSEED,IFRST C DATA JSEED,IFRST/123456789,0/ C END Subroutine ranorND(b,nd,ilosc_atomow2,ii1,seed1) double precision b C generates a random vector on a unit sphere onC an ND-dimensional space. C algorithm 40 page 410. C the algorithm assumes that the formula for bz is valid, i.e C bz=(1.0-2.0*ransq), C and also that, ransq can be computed as as the sum of the independent random C numbers to the second power (expression 3). dimension ran(nd-1), b(nd) integer seed,seed1 c write(*,*) "In Subroutine ranorND(b,nd), nd=", nd c write(*,*) "ilosc_atomow2=", ilosc_atomow2 c do j=1,ilosc_atomow2 ransq=2 seed = seed1*ii1 call srand(seed) do while (ransq.ge.1) c write(*,*) ransq rsq=0.0 do i=1,nd-1 ran(i)=1.0-2.0*rand(0) c write(*,*) i, ran(i) rsq=ran(i)*ran(i)+rsq ! (3) c write(*,*) i, rsq enddo ransq=rsq c write(*,*) j, ransq enddo ranh=2.0*sqrt(1.0-ransq) c write(*,*) j, ranh, ransq do i=1,nd-1 b(i)=ran(i)*ranh c write(*,*) ran(i), ranh enddo b(nd)=(1.0-2.0*ransq) c write(*,*) b(1,j), b(2,j), b(3,j), b(4,j) c enddo return end subroutine rot_trans(xt2,t,q,xf,ilosc_atomow2) C x coordinates, t translation vector, q quaternion double precision xt2(90000,3),xf(90000,3) double precision t(3),q(4) integer i write(*,*) "Rotation vector:", q(1), q(2), q(3), q(4) r11= q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4) r22= q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4) r33= q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4) r12= 2*(q(2)*q(3)-q(1)*q(4)) r21= 2*(q(2)*q(3)+q(1)*q(4)) r13= 2*(q(2)*q(4)+q(1)*q(3)) r31= 2*(q(2)*q(4)-q(1)*q(3)) r23= 2*(q(3)*q(4)-q(1)*q(2)) r32= 2*(q(3)*q(4)+q(1)*q(2)) c write(*,*) r11, r22, r33, r12, r21, r13, r31, r23, r32 do i=1,ilosc_atomow2 xf(i,1)= r11*xt2(i,1)+r12*xt2(i,2)+r13*xt2(i,3) + t(1) xf(i,2)= r21*xt2(i,1)+r22*xt2(i,2)+r23*xt2(i,3) + t(2) xf(i,3)= r31*xt2(i,1)+r32*xt2(i,2)+r33*xt2(i,3) + t(3) c write(*,*) i, xt2(i,1), xt2(i,2), xt2(i,3) c write(*,*) i, xf(i,1), xf(i,2), xf(i,3) enddo return end