examples/UNRES-Dock/generator_v13b.f

   1       program generator
   2       implicit none
   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
  19       character*54 linia
  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)
  22       real*8 xf(90000,3)
  23       character*1 chain
  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)
  28
  29       korekta1 = 1
  30       korekta2 = 1
  31
  32 c      write(*,*)"Sposob uruchamiania:"
  33 c      write(*,*)"nazwa_programu  nazwa_pliku_inp_pdb1
  34 c     &  nazwa_pliku_inp_pdb2  nazwa_pliku_out"
  35 c      write(*,*)""
  36
  37       nnn = iargc()
  38       if (nnn.lt.3) then
  39       write(*,*)"program  input_pdb1  input_pdb2
  40      &  [0-1]     0 for protein, 1 for peptide"
  41       stop
  42       endif
  43
  44       call getarg(1,plik1)
  45       call getarg(2,plik2)
  46       call getarg(3,plik3)
  47       open(1,file=plik1,status="old")
  48       open(2,file=plik2,status="old")
  49       open(4,file='temp1')
  50       open(5,file='temp2')
  51       write(*,*) "All files read"
  52
  53       read(plik3(1:1),'(i1)') peptide
  54
  55       if (peptide.ne.0.and.peptide.ne.1) then
  56         write(*,*) "Third parameter must be 0 or 1,
  57      &              where 0 is for a protein-protein docking,
  58      &              and 1 is for protein-peptide docking."
  59         STOP
  60       endif
  61
  62       write(*,*) "Option", peptide, "selected
  63      &            (0 for protein, 1 for peptide)"
  64
  65
  66 c MOLEKULA 1
  67       do while (.true.)
  68         read(1,'(a54)',end=22) linia
  69         if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then
  70         write(4,'(a54)') linia
  71         endif
  72         if (linia(1:3).eq.'TER') then
  73 c        write(2,'(a54)') linia
  74          goto 22
  75         endif
  76       enddo
  77
  78 22    continue
  79       close(4)
  80       write(*,*) "Molecule 1 read"
  81
  82
  83 c MOLEKULA 2
  84       do while (.true.)
  85         read(2,'(a54)',end=23) linia
  86         if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then
  87         write(5,'(a54)') linia
  88         endif
  89         if (linia(1:3).eq.'TER') then
  90 c        write(2,'(a54)') linia
  91          goto 23
  92         endif
  93       enddo
  94 23    continue
  95       close(5)
  96       write(*,*) "Molecule 2 read"
  97
  98       open(4,file='temp1')
  99       open(5,file='temp2')
 100
 101       minx1=999.99
 102       miny1=999.99
 103       minz1=999.99
 104
 105       maxx1=-999.99
 106       maxy1=-999.99
 107       maxz1=-999.99
 108
 109
 110       minx2=999.99
 111       miny2=999.99
 112       minz2=999.99
 113
 114       maxx2=-999.99
 115       maxy2=-999.99
 116       maxz2=-999.99
 117
 118
 119
 120 c MOLEKULA 1
 121       i=1
 122       xcm1=0
 123       ycm1=0
 124       zcm1=0
 125       ile_reszt_mol1=0
 126
 127       do while (.true.)
 128           read(4,800,end=111)
 129      &  atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),chain,numerr1(i),
 130      &  x1(i),y1(i),z1(i)
 131
 132
 133 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 134
 135          if (x1(i).lt.minx1) minx1=x1(i)
 136          if (y1(i).lt.miny1) miny1=y1(i)
 137          if (z1(i).lt.minz1) minz1=z1(i)
 138
 139          if (x1(i).gt.maxx1) maxx1=x1(i)
 140          if (y1(i).gt.maxy1) maxy1=y1(i)
 141          if (z1(i).gt.maxz1) maxz1=z1(i)
 142
 143 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 144
 145 c         write(*,*) nazwaa1(i)
 146          if (nazwaa1(i).eq.' CA ') then
 147              xcm1 = xcm1 + x1(i)
 148              ycm1 = ycm1 + y1(i)
 149              zcm1 = zcm1 + z1(i)
 150              ile_reszt_mol1 = ile_reszt_mol1 +1
 151          endif
 152
 153          if (atoma1(i).eq.'TER ') goto 111
 154          if (atoma1(i).eq.'END ') goto 111
 155          if (i.eq.1) write(3,'(a6)') "REMARK"
 156 c        write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),
 157 c     &               numerr1(i),x1(i),y1(i),z1(i)
 158
 159 c        write(*,*) i, x(i),y(i),z(i)
 160
 161         i=i+1
 162       enddo
 163 111   continue
 164
 165       distx1=maxx1-minx1
 166       disty1=maxy1-miny1
 167       distz1=maxz1-minz1
 168       ilosc_atomow1 = i-1
 169       xcm1 = xcm1 / (ile_reszt_mol1*1.0)
 170       ycm1 = ycm1 / (ile_reszt_mol1*1.0)
 171       zcm1 = zcm1 / (ile_reszt_mol1*1.0)
 172
 173       do i=1,ilosc_atomow1
 174          x1(i)=x1(i)-xcm1
 175          y1(i)=y1(i)-ycm1
 176          z1(i)=z1(i)-zcm1
 177       enddo
 178
 179       xcm1=0
 180       ycm1=0
 181       zcm1=0
 182
 183       do i=1,ilosc_atomow1
 184          if (nazwaa1(i).eq.' CA ') then
 185              xcm1 = xcm1 + x1(i)
 186              ycm1 = ycm1 + y1(i)
 187              zcm1 = zcm1 + z1(i)
 188          endif
 189       enddo
 190       xcm1 = xcm1 / (ile_reszt_mol1*1.0)
 191       ycm1 = ycm1 / (ile_reszt_mol1*1.0)
 192       zcm1 = zcm1 / (ile_reszt_mol1*1.0)
 193
 194       dist1cm=999.9
 195       dist1acm=0.0
 196       do i=1,ilosc_atomow1
 197 c        write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),
 198 c     &               numerr1(i),x1(i)-xcm1,y1(i)-ycm1,z1(i)-zcm1
 199
 200         dist1cmt = dsqrt(((xcm1-x1(i))**2)+
 201      &               ((ycm1-y1(i))**2)+
 202      &               ((zcm1-z1(i))**2))
 203
 204 c        write(*,*) dist1cmt, dist1cm, i, cm1
 205         if (dist1cm.gt.dist1cmt) then
 206            dist1cm=dist1cmt
 207            cm1=i
 208         endif
 209
 210       enddo
 211
 212       do i=1,ilosc_atomow1
 213       do j=i,ilosc_atomow1
 214                  dist1cmt = dsqrt(((x1(j)-x1(i))**2)+
 215      &               ((y1(j)-y1(i))**2)+
 216      &               ((z1(j)-z1(i))**2))
 217         if (dist1acm.lt.dist1cmt) then
 218            dist1acm=dist1cmt
 219         endif
 220       enddo
 221       enddo
 222
 223
 224 c      write(3,303) "TER"
 225 c      write(3,304) "ENDMDL"
 226
 227       write(*,*) "Molecule 1 calculated"
 228
 229 c MOLEKULA 2
 230       i=1
 231       xcm2=0
 232       ycm2=0
 233       zcm2=0
 234       ile_reszt_mol2=0
 235
 236       do while (.true.)
 237           read(5,800,end=112)
 238      &  atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),chain,numerr2(i),
 239      &  x2(i),y2(i),z2(i)
 240
 241 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 242
 243          if (x2(i).lt.minx2) minx2=x2(i)
 244          if (y2(i).lt.miny2) miny2=y2(i)
 245          if (z2(i).lt.minz2) minz2=z2(i)
 246
 247          if (x2(i).gt.maxx2) maxx2=x2(i)
 248          if (y2(i).gt.maxy2) maxy2=y2(i)
 249          if (z2(i).gt.maxz2) maxz2=z2(i)
 250
 251 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 252
 253 c         write(*,*) nazwaa1(i)
 254          if (nazwaa2(i).eq.' CA ') then
 255              xcm2 = xcm2 + x2(i)
 256              ycm2 = ycm2 + y2(i)
 257              zcm2 = zcm2 + z2(i)
 258              ile_reszt_mol2 = ile_reszt_mol2 + 1
 259          endif
 260
 261          if (atoma2(i).eq.'TER ') goto 111
 262          if (atoma2(i).eq.'END ') goto 111
 263 c         if (i.eq.1) write(3,'(a6)') "REMARK"
 264 c        write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),
 265 c     &               numerr2(i),x2(i),y2(i),z2(i)
 266 c        write(*,*) i, x(i),y(i),z(i)
 267
 268         i=i+1
 269       enddo
 270 112   continue
 271
 272       distx2=maxx2-minx2
 273       disty2=maxy2-miny2
 274       distz2=maxz2-minz2
 275       ilosc_atomow2 = i-1
 276
 277       xcm2 = xcm2 / (ile_reszt_mol2*1.0)
 278       ycm2 = ycm2 / (ile_reszt_mol2*1.0)
 279       zcm2 = zcm2 / (ile_reszt_mol2*1.0)
 280
 281       do i=1,ilosc_atomow2
 282          x2(i)=x2(i)-xcm2
 283          y2(i)=y2(i)-ycm2
 284          z2(i)=z2(i)-zcm2
 285 c        write(*,*) x2(i),y2(i),z2(i)
 286          xt2(i,1)=x2(i)
 287          xt2(i,2)=y2(i)
 288          xt2(i,3)=z2(i)
 289 c        write(*,*) xt2(i,1), xt2(i,2), xt2(i,3)
 290       enddo
 291
 292       xcm2=0
 293       ycm2=0
 294       zcm2=0
 295       do i=1,ilosc_atomow2
 296          if (nazwaa2(i).eq.' CA ') then
 297              xcm2 = xcm2 + x2(i)
 298              ycm2 = ycm2 + y2(i)
 299              zcm2 = zcm2 + z2(i)
 300          endif
 301       enddo
 302       xcm2 = xcm2 / (ile_reszt_mol2*1.0)
 303       ycm2 = ycm2 / (ile_reszt_mol2*1.0)
 304       zcm2 = zcm2 / (ile_reszt_mol2*1.0)
 305
 306
 307       dist2cm=999.9
 308       dist2acm=0.0
 309       do i=1,ilosc_atomow2
 310 c        write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),
 311 c     &               numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2
 312
 313         dist2cmt = dsqrt(((xcm2-x2(i))**2)+
 314      &               ((ycm2-y2(i))**2)+
 315      &               ((zcm2-z2(i))**2))
 316
 317 c        write(*,*) dist2cmt, dist2cm, i, cm2
 318         if (dist2cm.gt.dist2cmt) then
 319            dist2cm=dist2cmt
 320            cm2=i
 321         endif
 322
 323       enddo
 324
 325
 326       do i=1,ilosc_atomow2
 327       do j=i,ilosc_atomow2
 328                  dist2cmt = dsqrt(((x2(j)-x2(i))**2)+
 329      &               ((y2(j)-y2(i))**2)+
 330      &               ((z2(j)-z2(i))**2))
 331         if (dist2acm.lt.dist2cmt) then
 332            dist2acm=dist2cmt
 333         endif
 334       enddo
 335       enddo
 336
 337
 338 c      write(3,303) "TER"
 339
 340       glys1=0
 341       glys2=0
 342       glye1=0
 343       glye2=0
 344       r1s=nazwar1(1)
 345       r2s=nazwar2(1)
 346       r1e=nazwar1(ile_reszt_mol1)
 347       r2e=nazwar2(ile_reszt_mol2)
 348
 349       if (r1s.eq."GLY") then
 350         glys1=1
 351         write(*,*) "Chain 1 is starting from Gly ", r1s, glys1
 352       endif
 353
 354       if (r2s.eq."GLY") then
 355         glys2=1
 356         write(*,*) "Chain 2 is starting from Gly ", r2s, glys2
 357       endif
 358
 359       if (r1e.eq."GLY") then
 360         glye1=1
 361         write(*,*) "Chain 1 is ending at Gly ", r1e, glye1
 362       endif
 363
 364
 365       if (r2e.eq."GLY") then
 366         glye2=1
 367         write(*,*) "Chain 2 is ending from Gly ", r2e, glye2
 368       endif
 369
 370
 371 c checking for Gly
 372       ilee=4
 373       if(glys1.eq.1) ilee=ilee-1
 374       if(glys2.eq.1) ilee=ilee-1
 375       if(glye1.eq.1) ilee=ilee-1
 376       if(glye2.eq.1) ilee=ilee-1
 377       write(*,*) "How many dummy atoms:", ilee
 378       write(*,*) ""
 379
 380
 381       write(*,*) "Molecule 2 calculated"
 382       write(*,*) "Number of atoms:", ilosc_atomow1, ilosc_atomow2
 383
 384       write(*,*) "Min and max coordinates in the systems"
 385       write(*,300) minx1, maxx1, miny1, maxy1, minz1, maxz1
 386       write(*,300) minx2, maxx2, miny2, maxy2, minz2, maxz2
 387
 388       write(*,*) "Distances in axes x, y, z"
 389       write(*,301) distx1, disty1, distz1
 390       write(*,301) distx2, disty2, distz2
 391
 392       boxx = (distx1 + distx2)*1.5 + 20.0
 393       boxy = (disty1 + disty2)*1.5 + 20.0
 394       boxz = (distz1 + distz2)*1.5 + 20.0
 395
 396 c      write(*,*) "Boxsize [x, y, z]"
 397 c      write(*,301) boxx, boxy, boxz
 398
 399       write(*,*) "Maximum distances:"
 400       write(*,*) dist1acm, dist2acm
 401
 402       write(*,*) "Boxsize (suggested)"
 403       write(*,*) (dist1acm+dist2acm)*1.2+20
 404       open(61,file='boxsize.txt')
 405       write(61,302) ((dist1acm+dist2acm)*1.2+20)
 406       close (61)
 407
 408       write(*,*) "Center of masses:"
 409       write(*,301) xcm1, ycm1, zcm1
 410       write(*,301) xcm2, ycm2, zcm2
 411
 412       write(*,*) "Center of masses atoms:"
 413       write(*,800) atoma1(cm1),numera1(cm1),nazwaa1(cm1),nazwar1(cm1),
 414      &             " ",numerr1(cm1),x1(cm1),y1(cm1),z1(cm1)
 415       write(*,800) atoma2(cm2),numera2(cm2),nazwaa2(cm2),nazwar2(cm2),
 416      &             " ",numerr2(cm2),x2(cm2),y2(cm2),z2(cm2)
 417
 418
 419       write(*,*) "Number of amino-acid residues in two proteins:"
 420       write(*,*) ile_reszt_mol1, ile_reszt_mol2
 421
 422
 423 c ---------------------------------------------------------------------
 424 c KONIEC WCZYTYWANIA
 425 c ---------------------------------------------------------------------
 426
 427       korekta1 = numerr1(1)
 428       write(*,*) "Original number of the first residue in
 429      &            the first protein:"
 430       write(*,*) numerr1(1)
 431
 432       write(*,*) "Correction of the first residue in the first protein:"
 433       write(*,*) korekta1
 434
 435
 436 c       wymiar=(distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2)
 437 c       wymiar=wymiar/12.0
 438 c       wymiar=10.0
 439
 440 c      write(*,*) "Minimum required distance between centers of mass:"
 441 c      write(*,*) wymiar
 442
 443       boxx = (7+distx1 + dist2acm)/2.0
 444       boxy = (7+disty1 + dist2acm)/2.0
 445       boxz = (7+distz1 + dist2acm)/2.0
 446       write(*,*) "Rotation space (+/-):", boxx, boxy, boxz
 447
 448
 449       do ii1=1,24
 450
 451       proba = 1
 452       proba0x = 1
 453       proba0y = 1
 454       proba0z = 1
 455       proba1 = 1
 456       proba2 = 1
 457       proba3 = 1
 458       conf = 1
 459
 460
 461        if (ii1.lt.10) then
 462          write (plik3, "(A6,I1,a4)") "model0", ii1, ".pdb"
 463        else
 464          write (plik3, "(A5,I2,a4)") "model", ii1, ".pdb"
 465        endif
 466
 467        write(*,*) "Conformation will be saved to file:", plik3
 468
 469 c      plik3a = plik3//ilosc
 470       open(3,file=plik3)
 471
 472
 473       do i=1,ilosc_atomow1
 474         write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),"A",
 475      &               (numerr1(i)-korekta1+1),
 476      &               x1(i),y1(i),z1(i)
 477       enddo
 478       write(3,303) "TER"
 479 c      write(3,304) "ENDMDL"
 480
 481       if (peptide.eq.1) then
 482       wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/4
 483       endif
 484
 485       if (peptide.eq.0) then
 486       wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/7
 487       endif
 488
 489       if (ii1.eq.1) then
 490       write(*,*) "Initial minimum distance between centers of mass:",
 491      & wymiar
 492       endif
 493
 494 c ROTATIONS OF THE SECOND MOLECULE
 495       nd=4
 496       t(1)=0.0
 497       t(2)=0.0
 498       t(3)=0.0
 499       call ranorND(b,nd,ilosc_atomow2,ii1)
 500 c      do i=1,4
 501 c        write(*,*) b(i)
 502 c      enddo
 503       call rot_trans(xt2,t,b,xf,ilosc_atomow2)
 504       do i=1,ilosc_atomow2
 505         x2(i)=xf(i,1)
 506         y2(i)=xf(i,2)
 507         z2(i)=xf(i,3)
 508       enddo
 509 c END OF ROTATIONS
 510
 511
 512        do while (.true.)
 513         if (conf.gt.1) goto 133
 514
 515         seed = time()*proba
 516         call srand(seed)
 517         r = rand(0)
 518 c        write(*,*)rand(0),rand(0),rand(0)
 519 c        write(*,*)rand(0),rand(0),rand(0)
 520
 521         cm2xr = (rand(0)*2*boxx)-(boxx)
 522         cm2yr = (rand(0)*2*boxy)-(boxy)
 523         cm2zr = (rand(0)*2*boxz)-(boxz)
 524 c        write(*,*) cm2xr, cm2yr, cm2zr
 525 c        write(*,300)rand(0),rand(0),rand(0)
 526
 527 c        write(*,*) ABS(cm2xr+xcm1), (dist1acm+dist2acm+3)
 528
 529 c        if (ABS(cm2xr+xcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
 530 c           proba0x=proba0x+1
 531 c           goto 200
 532 c        endif
 533
 534 c        if (ABS(cm2yr+ycm1).gt.((dist1acm+dist2acm+5)*0.8)) then
 535 c           proba0y=proba0y+1
 536 c           goto 200
 537 c        endif
 538
 539 c        if (ABS(cm2zr+zcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
 540 c           proba0z=proba0z+1
 541 c           goto 200
 542 c        endif
 543
 544
 545       do i=1,ilosc_atomow2
 546         x2r(i)=x2(i)-cm2xr
 547         y2r(i)=y2(i)-cm2yr
 548         z2r(i)=z2(i)-cm2zr
 549       enddo
 550
 551
 552 c SPRAWDZANIE CZY NIE JEST ZBYT PODOBNE DO JUZ WYLOSOWANEGO
 553       cmxx(ii1)=cm2xr
 554       cmyy(ii1)=cm2yr
 555       cmzz(ii1)=cm2zr
 556
 557       if (ii1.gt.1) then
 558
 559        do ir=1,ii1-1
 560          dd=sqrt((cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+
 561      &          (cm2zr-cmzz(ir))**2)
 562
 563 c         write(*,*) "Center of mass is", dd, "A away from", ir
 564 c         write(*,*) "Minimum distance required:", wymiar, "A"
 565
 566 c           write(*,*) ii1, ir, proba3, dd, wymiar
 567
 568          if (dd.lt.wymiar) then
 569            proba3=proba3+1
 570          if (proba3.gt.1.and.peptide.eq.1) wymiar=wymiar*0.99992
 571          if (proba3.gt.1.and.peptide.eq.0) wymiar=wymiar*0.99992
 572 c           write(*,*) proba3, wymiar
 573            goto 200
 574          endif
 575
 576        enddo
 577       endif
 578
 579
 580
 581
 582 c SPRAWDZANIE ODLEGLOSCI
 583         accept=0
 584       do i=1,ilosc_atomow2
 585 c        write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),
 586 c     &               numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2
 587
 588         do j=1,ilosc_atomow1
 589            dist12 = sqrt(((x2r(i)-x1(j))**2)+
 590      &               ((y2r(i)-y1(j))**2)+
 591      &               ((z2r(i)-z1(j))**2))
 592 c           write(*,*) dist12, cm2xr, cm2yr, cm2zr
 593            if (dist12.lt.2.5) then
 594                proba1=proba1+1
 595                goto 200
 596            endif
 597
 598            if (dist12.le.3.5) accept=1
 599 c           if (dist12.le.3.5) accept=1
 600 c           if (dist12.gt.3.5) proba2=proba2+1
 601         enddo
 602       enddo
 603
 604 c      write(*,*) "No overlaps", proba
 605       if (accept.eq.0) goto 200
 606 c      write(*,*) "No overlaps, contact", conf, proba
 607 c      write(3,305) "MODEL", conf
 608 c      conf = conf+1
 609
 610
 611
 612        write(*,*) "Conformation", ii1, "Is accepted after",
 613      &             proba0x+proba0y+proba0z,
 614      &             proba1, proba2, proba3, proba, "attempts. ",
 615      &             "Dist between CoM of previous=", wymiar
 616
 617 c Dodatkowe info
 618       write(*,*)"Extr",ABS(cm2xr+xcm1), ABS(cm2yr+ycm1), ABS(cm2zr+zcm1)
 619
 620
 621        conf=conf+1
 622 c       write (*,*)
 623
 624
 625
 626
 627
 628       korekta2 = numerr2(1)-1
 629 c      write(*,*) "Original number of the first residue in
 630 c     &            the second protein:"
 631 c      write(*,*) numerr2(1)
 632
 633 c      write(*,*)"Correction of the first residue in the second protein:"
 634 c      write(*,*) korekta2
 635       do i=1,ilosc_atomow2
 636         write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),"B",
 637      &               (numerr2(i)+ile_reszt_mol1-korekta2),
 638      &               x2r(i),y2r(i),z2r(i)
 639       enddo
 640         write(3,303) "TER"
 641 c        write(3,304) "ENDMDL"
 642 c        write(3,305) "MODEL", conf
 643         write(3,303) "END"
 644
 645  200  continue
 646         proba = proba + 1
 647       enddo
 648  133  continue
 649
 650       close (3)
 651
 652 c      call sleep(1)
 653
 654
 655        if (ii1.lt.10) then
 656          write (plik4, "(A9,I1,a4)") "template0", ii1, ".sco"
 657        else
 658          write (plik4, "(A8,I2,a4)") "template", ii1, ".sco"
 659        endif
 660
 661        open (7,file=plik4)
 662
 663
 664       if (peptide.eq.1) then
 665       do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee)
 666         if (glys1.eq.0) then
 667           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 668           if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
 669      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 670           if (ijl.ge.(ile_reszt_mol1+2))
 671      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 672         endif
 673
 674         if (glys1.eq.1) then
 675           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 676           if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
 677      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 678           if (ijl.ge.(ile_reszt_mol1+1))
 679      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 680         endif
 681       enddo
 682       endif
 683
 684
 685       if (peptide.eq.0) then
 686       do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee)
 687 c        write(*,*) "ijl"
 688 c no Gly
 689         if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.0) then
 690           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 691           if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
 692      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 693           if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3))
 694      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 695           if (ijl.gt.(ile_reszt_mol1+3).and.
 696      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1+3))
 697      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 698           if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+4))
 699      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 700         endif
 701
 702
 703 c chain 1 first Gly
 704         if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.0) then
 705           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 706           if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
 707      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 708           if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+2))
 709      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 710           if (ijl.gt.(ile_reszt_mol1+2).and.
 711      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2))
 712      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 713           if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3))
 714      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 715         endif
 716
 717 c chain 1 last Gly
 718         if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.0) then
 719           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 720           if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
 721      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 722           if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2))
 723      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 724           if (ijl.gt.(ile_reszt_mol1+2).and.
 725      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2))
 726      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 727           if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3))
 728      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 729         endif
 730
 731 c chain 1 first Gly and chain 1 last Gly
 732         if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.0) then
 733           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 734           if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
 735      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 736           if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1))
 737      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 738           if (ijl.gt.(ile_reszt_mol1+1).and.
 739      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1))
 740      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 741           if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2))
 742      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 743         endif
 744
 745 c chain 2 first Gly
 746         if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.1) then
 747           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 748           if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
 749      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 750           if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2))
 751      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 752           if (ijl.gt.(ile_reszt_mol1+2).and.
 753      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2))
 754      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 755           if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3))
 756      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 757         endif
 758
 759 c chain 1 first Gly chain 2 first Gly
 760         if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.1) then
 761           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 762           if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
 763      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 764           if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1))
 765      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 766           if (ijl.gt.(ile_reszt_mol1+1).and.
 767      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1))
 768      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 769           if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2))
 770      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 771         endif
 772
 773 c chain 1 last Gly chain 2 first Gly
 774         if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.1) then
 775           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 776           if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1)
 777      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 778 c          if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3))
 779 c     &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 780           if (ijl.gt.(ile_reszt_mol1+1).and.
 781      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1))
 782      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 783           if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2))
 784      &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 785         endif
 786
 787 c all Gly
 788         if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.1) then
 789           if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 790           if (ijl.le.ile_reszt_mol1.and.ijl.gt.1)
 791      &      write(7,*) ijl, "1.0 1.0 1.0 ", "1"
 792 c          if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3))
 793 c     &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 794           if (ijl.gt.(ile_reszt_mol1).and.
 795      &     ijl.le.(ile_reszt_mol2+ile_reszt_mol1))
 796      &      write(7,*) ijl, "1.0 1.0 1.0 ", "2"
 797 c          if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+1))
 798 c     &      write(7,*) ijl, "1.0 1.0 1.0 ", "0"
 799         endif
 800
 801
 802
 803       enddo
 804       endif
 805
 806
 807        close (7)
 808
 809 c     koniec petli po ilosci struktur
 810       enddo
 811
 812 c      do i=1,20
 813 c        write(*,*) i, cmxx(i), cmyy(i), cmzz(i)
 814 c         dd=(cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+(cm2zr-cmyy(ir))**2
 815 c      enddo
 816
 817 c      do i=1,20
 818 c      do j=1,20
 819 c         dd=sqrt((cmxx(i)-cmxx(j))**2+(cmyy(i)-cmyy(j))**2+
 820 c     &      (cmzz(i)-cmzz(j))**2)
 821 c         write(*,*) i,j,dd
 822 c      enddo
 823 c      enddo
 824
 825
 826   300 format(6f8.2)
 827   301 format(3f8.2)
 828   302 format(f8.2)
 829   800 format(a4,2X,i5,1X,a4,1X,a3,1X,a1,i4,4X,3f8.3)
 830   303 format(a3)
 831   304 format(a6)
 832   305 format(a5,1X,i8)
 833       end
 834
 835
 836
 837
 838       SUBROUTINE SRAND(ISEED)
 839 C
 840 C  This subroutine sets the integer seed to be used with the
 841 C  companion RAND function to the value of ISEED.  A flag is
 842 C  set to indicate that the sequence of pseudo-random numbers
 843 C  for the specified seed should start from the beginning.
 844 C
 845       COMMON /SEED/JSEED,IFRST
 846 C
 847       JSEED = ISEED
 848       IFRST = 0
 849 C
 850       RETURN
 851       END
 852       REAL FUNCTION RAND()
 853 C
 854 C  This function returns a pseudo-random number for each invocation.
 855 C  It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal
 856 C  standard number generator whose Pascal code appears in the article:
 857 C
 858 C     Park, Steven K. and Miller, Keith W., "Random Number Generators:
 859 C     Good Ones are Hard to Find", Communications of the ACM,
 860 C     October, 1988.
 861 C
 862       PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773,
 863      +           MOMDMP=2836)
 864 C
 865       COMMON  /SEED/JSEED,IFRST
 866       INTEGER HVLUE, LVLUE, TESTV, NEXTN
 867       SAVE    NEXTN
 868 C
 869       IF (IFRST .EQ. 0) THEN
 870         NEXTN = JSEED
 871         IFRST = 1
 872       ENDIF
 873 C
 874       HVLUE = NEXTN / MOBYMP
 875       LVLUE = MOD(NEXTN, MOBYMP)
 876       TESTV = MPLIER*LVLUE - MOMDMP*HVLUE
 877       IF (TESTV .GT. 0) THEN
 878         NEXTN = TESTV
 879       ELSE
 880         NEXTN = TESTV + MODLUS
 881       ENDIF
 882       RAND = REAL(NEXTN)/REAL(MODLUS)
 883 C
 884       RETURN
 885       END
 886       BLOCKDATA RANDBD
 887       COMMON /SEED/JSEED,IFRST
 888 C
 889       DATA JSEED,IFRST/123456789,0/
 890 C
 891       END
 892
 893       Subroutine ranorND(b,nd,ilosc_atomow2,ii1)
 894       double precision b
 895 C generates a random vector on a unit sphere onC an ND-dimensional space.
 896 C algorithm 40 page 410.
 897 C  the algorithm assumes that the formula for bz is valid, i.e
 898 C   bz=(1.0-2.0*ransq),
 899 C and also that, ransq can be computed as  as the sum of the independent random
 900 C numbers to the second power (expression 3).
 901       dimension ran(nd-1), b(nd)
 902       integer seed
 903 c      write(*,*) "In Subroutine ranorND(b,nd), nd=", nd
 904 c      write(*,*) "ilosc_atomow2=", ilosc_atomow2
 905 c      do j=1,ilosc_atomow2
 906         ransq=2
 907         seed = time()*ii1
 908         call srand(seed)
 909         do while (ransq.ge.1)
 910 c       write(*,*) ransq
 911            rsq=0.0
 912            do i=1,nd-1
 913               ran(i)=1.0-2.0*rand(0)
 914 c           write(*,*) i, ran(i)
 915               rsq=ran(i)*ran(i)+rsq          !  (3)
 916 c           write(*,*) i, rsq
 917            enddo
 918            ransq=rsq
 919 c           write(*,*) j, ransq
 920         enddo
 921         ranh=2.0*sqrt(1.0-ransq)
 922 c       write(*,*) j, ranh, ransq
 923         do i=1,nd-1
 924           b(i)=ran(i)*ranh
 925 c          write(*,*) ran(i), ranh
 926         enddo
 927         b(nd)=(1.0-2.0*ransq)
 928 c        write(*,*) b(1,j), b(2,j), b(3,j), b(4,j)
 929 c      enddo
 930       return
 931       end
 932
 933       subroutine rot_trans(xt2,t,q,xf,ilosc_atomow2)
 934 C x coordinates, t translation vector, q quaternion
 935       double precision xt2(90000,3),xf(90000,3)
 936       double precision t(3),q(4)
 937       integer i
 938       write(*,*) "Rotation vector:", q(1), q(2), q(3), q(4)
 939       r11= q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4)
 940       r22= q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4)
 941       r33= q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4)
 942       r12= 2*(q(2)*q(3)-q(1)*q(4))
 943       r21= 2*(q(2)*q(3)+q(1)*q(4))
 944       r13= 2*(q(2)*q(4)+q(1)*q(3))
 945       r31= 2*(q(2)*q(4)-q(1)*q(3))
 946       r23= 2*(q(3)*q(4)-q(1)*q(2))
 947       r32= 2*(q(3)*q(4)+q(1)*q(2))
 948
 949 c      write(*,*) r11, r22, r33, r12, r21, r13, r31, r23, r32
 950
 951       do i=1,ilosc_atomow2
 952         xf(i,1)=  r11*xt2(i,1)+r12*xt2(i,2)+r13*xt2(i,3) + t(1)
 953         xf(i,2)=  r21*xt2(i,1)+r22*xt2(i,2)+r23*xt2(i,3) + t(2)
 954         xf(i,3)=  r31*xt2(i,1)+r32*xt2(i,2)+r33*xt2(i,3) + t(3)
 955 c        write(*,*) i, xt2(i,1), xt2(i,2), xt2(i,3)
 956 c        write(*,*) i, xf(i,1), xf(i,2), xf(i,3)
 957       enddo
 958
 959       return
 960       end
 961