--- /dev/null
+ 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
+ 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
+ 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)
+ integer time_a(3)
+
+ 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.3) then
+ write(*,*)"program input_pdb1 input_pdb2
+ & [0-1] 0 for protein, 1 for peptide"
+ stop
+ endif
+
+ call getarg(1,plik1)
+ call getarg(2,plik2)
+ call getarg(3,plik3)
+ 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
+
+ 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 = (7+distx1 + dist2acm)/2.0
+ boxy = (7+disty1 + dist2acm)/2.0
+ boxz = (7+distz1 + dist2acm)/2.0
+ write(*,*) "Rotation space (+/-):", boxx, boxy, boxz
+
+
+ do ii1=1,24
+
+ 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))/7
+ 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)
+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
+
+#ifdef GFORTRAN
+ seed = time()*proba
+#else
+ call itime(time_a)
+ seed = (60*time_a(2)+time_a(3))*proba
+#endif
+ call srand(seed)
+ r = rand(0)
+c write(*,*)rand(0),rand(0),rand(0)
+c write(*,*)rand(0),rand(0),rand(0)
+
+ cm2xr = (rand(0)*2*boxx)-(boxx)
+ cm2yr = (rand(0)*2*boxy)-(boxy)
+ cm2zr = (rand(0)*2*boxz)-(boxz)
+c write(*,*) cm2xr, cm2yr, cm2zr
+c write(*,300)rand(0),rand(0),rand(0)
+
+c write(*,*) ABS(cm2xr+xcm1), (dist1acm+dist2acm+3)
+
+c if (ABS(cm2xr+xcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
+c proba0x=proba0x+1
+c goto 200
+c endif
+
+c if (ABS(cm2yr+ycm1).gt.((dist1acm+dist2acm+5)*0.8)) then
+c proba0y=proba0y+1
+c goto 200
+c endif
+
+c if (ABS(cm2zr+zcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
+c proba0z=proba0z+1
+c goto 200
+c endif
+
+
+ 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)
+ 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
+c write(*,*) "In Subroutine ranorND(b,nd), nd=", nd
+c write(*,*) "ilosc_atomow2=", ilosc_atomow2
+c do j=1,ilosc_atomow2
+ ransq=2
+ seed = time()*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
+
projects / unres.git / commitdiff