source/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       integer time_a(3)
  29
  30       korekta1 = 1
  31       korekta2 = 1
  32
  33 c      write(*,*)"Sposob uruchamiania:"
  34 c      write(*,*)"nazwa_programu  nazwa_pliku_inp_pdb1
  35 c     &  nazwa_pliku_inp_pdb2  nazwa_pliku_out"
  36 c      write(*,*)""
  37
  38       nnn = iargc()
  39       if (nnn.lt.3) then
  40       write(*,*)"program  input_pdb1  input_pdb2
  41      &  [0-1]     0 for protein, 1 for peptide"
  42       stop
  43       endif
  44
  45       call getarg(1,plik1)
  46       call getarg(2,plik2)
  47       call getarg(3,plik3)
  48       open(1,file=plik1,status="old")
  49       open(2,file=plik2,status="old")
  50       open(4,file='temp1')
  51       open(5,file='temp2')
  52       write(*,*) "All files read"
  53
  54       read(plik3(1:1),'(i1)') peptide
  55
  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."
  60         STOP
  61       endif
  62
  63       write(*,*) "Option", peptide, "selected
  64      &            (0 for protein, 1 for peptide)"
  65
  66
  67 c MOLEKULA 1
  68       do while (.true.)
  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
  72         endif
  73         if (linia(1:3).eq.'TER') then
  74 c        write(2,'(a54)') linia
  75          goto 22
  76         endif
  77       enddo
  78
  79 22    continue
  80       close(4)
  81       write(*,*) "Molecule 1 read"
  82
  83
  84 c MOLEKULA 2
  85       do while (.true.)
  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
  89         endif
  90         if (linia(1:3).eq.'TER') then
  91 c        write(2,'(a54)') linia
  92          goto 23
  93         endif
  94       enddo
  95 23    continue
  96       close(5)
  97       write(*,*) "Molecule 2 read"
  98
  99       open(4,file='temp1')
 100       open(5,file='temp2')
 101
 102       minx1=999.99
 103       miny1=999.99
 104       minz1=999.99
 105
 106       maxx1=-999.99
 107       maxy1=-999.99
 108       maxz1=-999.99
 109
 110
 111       minx2=999.99
 112       miny2=999.99
 113       minz2=999.99
 114
 115       maxx2=-999.99
 116       maxy2=-999.99
 117       maxz2=-999.99
 118
 119
 120
 121 c MOLEKULA 1
 122       i=1
 123       xcm1=0
 124       ycm1=0
 125       zcm1=0
 126       ile_reszt_mol1=0
 127
 128       do while (.true.)
 129           read(4,800,end=111)
 130      &  atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),chain,numerr1(i),
 131      &  x1(i),y1(i),z1(i)
 132
 133
 134 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 135
 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)
 139
 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)
 143
 144 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 145
 146 c         write(*,*) nazwaa1(i)
 147          if (nazwaa1(i).eq.' CA ') then
 148              xcm1 = xcm1 + x1(i)
 149              ycm1 = ycm1 + y1(i)
 150              zcm1 = zcm1 + z1(i)
 151              ile_reszt_mol1 = ile_reszt_mol1 +1
 152          endif
 153
 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)
 159
 160 c        write(*,*) i, x(i),y(i),z(i)
 161
 162         i=i+1
 163       enddo
 164 111   continue
 165
 166       distx1=maxx1-minx1
 167       disty1=maxy1-miny1
 168       distz1=maxz1-minz1
 169       ilosc_atomow1 = i-1
 170       xcm1 = xcm1 / (ile_reszt_mol1*1.0)
 171       ycm1 = ycm1 / (ile_reszt_mol1*1.0)
 172       zcm1 = zcm1 / (ile_reszt_mol1*1.0)
 173
 174       do i=1,ilosc_atomow1
 175          x1(i)=x1(i)-xcm1
 176          y1(i)=y1(i)-ycm1
 177          z1(i)=z1(i)-zcm1
 178       enddo
 179
 180       xcm1=0
 181       ycm1=0
 182       zcm1=0
 183
 184       do i=1,ilosc_atomow1
 185          if (nazwaa1(i).eq.' CA ') then
 186              xcm1 = xcm1 + x1(i)
 187              ycm1 = ycm1 + y1(i)
 188              zcm1 = zcm1 + z1(i)
 189          endif
 190       enddo
 191       xcm1 = xcm1 / (ile_reszt_mol1*1.0)
 192       ycm1 = ycm1 / (ile_reszt_mol1*1.0)
 193       zcm1 = zcm1 / (ile_reszt_mol1*1.0)
 194
 195       dist1cm=999.9
 196       dist1acm=0.0
 197       do i=1,ilosc_atomow1
 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
 200
 201         dist1cmt = dsqrt(((xcm1-x1(i))**2)+
 202      &               ((ycm1-y1(i))**2)+
 203      &               ((zcm1-z1(i))**2))
 204
 205 c        write(*,*) dist1cmt, dist1cm, i, cm1
 206         if (dist1cm.gt.dist1cmt) then
 207            dist1cm=dist1cmt
 208            cm1=i
 209         endif
 210
 211       enddo
 212
 213       do i=1,ilosc_atomow1
 214       do j=i,ilosc_atomow1
 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
 219            dist1acm=dist1cmt
 220         endif
 221       enddo
 222       enddo
 223
 224
 225 c      write(3,303) "TER"
 226 c      write(3,304) "ENDMDL"
 227
 228       write(*,*) "Molecule 1 calculated"
 229
 230 c MOLEKULA 2
 231       i=1
 232       xcm2=0
 233       ycm2=0
 234       zcm2=0
 235       ile_reszt_mol2=0
 236
 237       do while (.true.)
 238           read(5,800,end=112)
 239      &  atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),chain,numerr2(i),
 240      &  x2(i),y2(i),z2(i)
 241
 242 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 243
 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)
 247
 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)
 251
 252 c         write(*,300) minx, maxx, miny, maxy, minz, maxz
 253
 254 c         write(*,*) nazwaa1(i)
 255          if (nazwaa2(i).eq.' CA ') then
 256              xcm2 = xcm2 + x2(i)
 257              ycm2 = ycm2 + y2(i)
 258              zcm2 = zcm2 + z2(i)
 259              ile_reszt_mol2 = ile_reszt_mol2 + 1
 260          endif
 261
 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)
 268
 269         i=i+1
 270       enddo
 271 112   continue
 272
 273       distx2=maxx2-minx2
 274       disty2=maxy2-miny2
 275       distz2=maxz2-minz2
 276       ilosc_atomow2 = i-1
 277
 278       xcm2 = xcm2 / (ile_reszt_mol2*1.0)
 279       ycm2 = ycm2 / (ile_reszt_mol2*1.0)
 280       zcm2 = zcm2 / (ile_reszt_mol2*1.0)
 281
 282       do i=1,ilosc_atomow2
 283          x2(i)=x2(i)-xcm2
 284          y2(i)=y2(i)-ycm2
 285          z2(i)=z2(i)-zcm2
 286 c        write(*,*) x2(i),y2(i),z2(i)
 287          xt2(i,1)=x2(i)
 288          xt2(i,2)=y2(i)
 289          xt2(i,3)=z2(i)
 290 c        write(*,*) xt2(i,1), xt2(i,2), xt2(i,3)
 291       enddo
 292
 293       xcm2=0
 294       ycm2=0
 295       zcm2=0
 296       do i=1,ilosc_atomow2
 297          if (nazwaa2(i).eq.' CA ') then
 298              xcm2 = xcm2 + x2(i)
 299              ycm2 = ycm2 + y2(i)
 300              zcm2 = zcm2 + z2(i)
 301          endif
 302       enddo
 303       xcm2 = xcm2 / (ile_reszt_mol2*1.0)
 304       ycm2 = ycm2 / (ile_reszt_mol2*1.0)
 305       zcm2 = zcm2 / (ile_reszt_mol2*1.0)
 306
 307
 308       dist2cm=999.9
 309       dist2acm=0.0
 310       do i=1,ilosc_atomow2
 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
 313
 314         dist2cmt = dsqrt(((xcm2-x2(i))**2)+
 315      &               ((ycm2-y2(i))**2)+
 316      &               ((zcm2-z2(i))**2))
 317
 318 c        write(*,*) dist2cmt, dist2cm, i, cm2
 319         if (dist2cm.gt.dist2cmt) then
 320            dist2cm=dist2cmt
 321            cm2=i
 322         endif
 323
 324       enddo
 325
 326
 327       do i=1,ilosc_atomow2
 328       do j=i,ilosc_atomow2
 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
 333            dist2acm=dist2cmt
 334         endif
 335       enddo
 336       enddo
 337
 338
 339 c      write(3,303) "TER"
 340
 341       glys1=0
 342       glys2=0
 343       glye1=0
 344       glye2=0
 345       r1s=nazwar1(1)
 346       r2s=nazwar2(1)
 347       r1e=nazwar1(ile_reszt_mol1)
 348       r2e=nazwar2(ile_reszt_mol2)
 349
 350       if (r1s.eq."GLY") then
 351         glys1=1
 352         write(*,*) "Chain 1 is starting from Gly ", r1s, glys1
 353       endif
 354
 355       if (r2s.eq."GLY") then
 356         glys2=1
 357         write(*,*) "Chain 2 is starting from Gly ", r2s, glys2
 358       endif
 359
 360       if (r1e.eq."GLY") then
 361         glye1=1
 362         write(*,*) "Chain 1 is ending at Gly ", r1e, glye1
 363       endif
 364
 365
 366       if (r2e.eq."GLY") then
 367         glye2=1
 368         write(*,*) "Chain 2 is ending from Gly ", r2e, glye2
 369       endif
 370
 371
 372 c checking for Gly
 373       ilee=4
 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
 379       write(*,*) ""
 380
 381
 382       write(*,*) "Molecule 2 calculated"
 383       write(*,*) "Number of atoms:", ilosc_atomow1, ilosc_atomow2
 384
 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
 388
 389       write(*,*) "Distances in axes x, y, z"
 390       write(*,301) distx1, disty1, distz1
 391       write(*,301) distx2, disty2, distz2
 392
 393       boxx = (distx1 + distx2)*1.5 + 20.0
 394       boxy = (disty1 + disty2)*1.5 + 20.0
 395       boxz = (distz1 + distz2)*1.5 + 20.0
 396
 397 c      write(*,*) "Boxsize [x, y, z]"
 398 c      write(*,301) boxx, boxy, boxz
 399
 400       write(*,*) "Maximum distances:"
 401       write(*,*) dist1acm, dist2acm
 402
 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)
 407       close (61)
 408
 409       write(*,*) "Center of masses:"
 410       write(*,301) xcm1, ycm1, zcm1
 411       write(*,301) xcm2, ycm2, zcm2
 412
 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)
 418
 419
 420       write(*,*) "Number of amino-acid residues in two proteins:"
 421       write(*,*) ile_reszt_mol1, ile_reszt_mol2
 422
 423
 424 c ---------------------------------------------------------------------
 425 c KONIEC WCZYTYWANIA
 426 c ---------------------------------------------------------------------
 427
 428       korekta1 = numerr1(1)
 429       write(*,*) "Original number of the first residue in
 430      &            the first protein:"
 431       write(*,*) numerr1(1)
 432
 433       write(*,*) "Correction of the first residue in the first protein:"
 434       write(*,*) korekta1
 435
 436
 437 c       wymiar=(distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2)
 438 c       wymiar=wymiar/12.0
 439 c       wymiar=10.0
 440
 441 c      write(*,*) "Minimum required distance between centers of mass:"
 442 c      write(*,*) wymiar
 443
 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
 448
 449
 450       do ii1=1,24
 451
 452       proba = 1
 453       proba0x = 1
 454       proba0y = 1
 455       proba0z = 1
 456       proba1 = 1
 457       proba2 = 1
 458       proba3 = 1
 459       conf = 1
 460
 461
 462        if (ii1.lt.10) then
 463          write (plik3, "(A6,I1,a4)") "model0", ii1, ".pdb"
 464        else
 465          write (plik3, "(A5,I2,a4)") "model", ii1, ".pdb"
 466        endif
 467
 468        write(*,*) "Conformation will be saved to file:", plik3
 469
 470 c      plik3a = plik3//ilosc
 471       open(3,file=plik3)
 472
 473
 474       do i=1,ilosc_atomow1
 475         write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),"A",
 476      &               (numerr1(i)-korekta1+1),
 477      &               x1(i),y1(i),z1(i)
 478       enddo
 479       write(3,303) "TER"
 480 c      write(3,304) "ENDMDL"
 481
 482       if (peptide.eq.1) then
 483       wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/4
 484       endif
 485
 486       if (peptide.eq.0) then
 487       wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/7
 488       endif
 489
 490       if (ii1.eq.1) then
 491       write(*,*) "Initial minimum distance between centers of mass:",
 492      & wymiar
 493       endif
 494
 495 c ROTATIONS OF THE SECOND MOLECULE
 496       nd=4
 497       t(1)=0.0
 498       t(2)=0.0
 499       t(3)=0.0
 500       call ranorND(b,nd,ilosc_atomow2,ii1)
 501 c      do i=1,4
 502 c        write(*,*) b(i)
 503 c      enddo
 504       call rot_trans(xt2,t,b,xf,ilosc_atomow2)
 505       do i=1,ilosc_atomow2
 506         x2(i)=xf(i,1)
 507         y2(i)=xf(i,2)
 508         z2(i)=xf(i,3)
 509       enddo
 510 c END OF ROTATIONS
 511
 512
 513        do while (.true.)
 514         if (conf.gt.1) goto 133
 515
 516 #ifdef GFORTRAN
 517         seed = time()*proba
 518 #else
 519         call  itime(time_a)
 520         seed = (60*time_a(2)+time_a(3))*proba
 521 #endif
 522         call srand(seed)
 523         r = rand(0)
 524 c        write(*,*)rand(0),rand(0),rand(0)
 525 c        write(*,*)rand(0),rand(0),rand(0)
 526
 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)
 532
 533 c        write(*,*) ABS(cm2xr+xcm1), (dist1acm+dist2acm+3)
 534
 535 c        if (ABS(cm2xr+xcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
 536 c           proba0x=proba0x+1
 537 c           goto 200
 538 c        endif
 539
 540 c        if (ABS(cm2yr+ycm1).gt.((dist1acm+dist2acm+5)*0.8)) then
 541 c           proba0y=proba0y+1
 542 c           goto 200
 543 c        endif
 544
 545 c        if (ABS(cm2zr+zcm1).gt.((dist1acm+dist2acm+5)*0.8)) then
 546 c           proba0z=proba0z+1
 547 c           goto 200
 548 c        endif
 549
 550
 551       do i=1,ilosc_atomow2
 552         x2r(i)=x2(i)-cm2xr
 553         y2r(i)=y2(i)-cm2yr
 554         z2r(i)=z2(i)-cm2zr
 555       enddo
 556
 557
 558 c SPRAWDZANIE CZY NIE JEST ZBYT PODOBNE DO JUZ WYLOSOWANEGO
 559       cmxx(ii1)=cm2xr
 560       cmyy(ii1)=cm2yr
 561       cmzz(ii1)=cm2zr
 562
 563       if (ii1.gt.1) then
 564
 565        do ir=1,ii1-1
 566          dd=sqrt((cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+
 567      &          (cm2zr-cmzz(ir))**2)
 568
 569 c         write(*,*) "Center of mass is", dd, "A away from", ir
 570 c         write(*,*) "Minimum distance required:", wymiar, "A"
 571
 572 c           write(*,*) ii1, ir, proba3, dd, wymiar
 573
 574          if (dd.lt.wymiar) then
 575            proba3=proba3+1
 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
 579            goto 200
 580          endif
 581
 582        enddo
 583       endif
 584
 585
 586
 587
 588 c SPRAWDZANIE ODLEGLOSCI
 589         accept=0
 590       do i=1,ilosc_atomow2
 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
 593
 594         do j=1,ilosc_atomow1
 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
 600                proba1=proba1+1
 601                goto 200
 602            endif
 603
 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
 607         enddo
 608       enddo
 609
 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
 614 c      conf = conf+1
 615
 616
 617
 618        write(*,*) "Conformation", ii1, "Is accepted after",
 619      &             proba0x+proba0y+proba0z,
 620      &             proba1, proba2, proba3, proba, "attempts. ",
 621      &             "Dist between CoM of previous=", wymiar
 622
 623 c Dodatkowe info
 624       write(*,*)"Extr",ABS(cm2xr+xcm1), ABS(cm2yr+ycm1), ABS(cm2zr+zcm1)
 625
 626
 627        conf=conf+1
 628 c       write (*,*)
 629
 630
 631
 632
 633
 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)
 638
 639 c      write(*,*)"Correction of the first residue in the second protein:"
 640 c      write(*,*) korekta2
 641       do i=1,ilosc_atomow2
 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)
 645       enddo
 646         write(3,303) "TER"
 647 c        write(3,304) "ENDMDL"
 648 c        write(3,305) "MODEL", conf
 649         write(3,303) "END"
 650
 651  200  continue
 652         proba = proba + 1
 653       enddo
 654  133  continue
 655
 656       close (3)
 657
 658 c      call sleep(1)
 659
 660
 661        if (ii1.lt.10) then
 662          write (plik4, "(A9,I1,a4)") "template0", ii1, ".sco"
 663        else
 664          write (plik4, "(A8,I2,a4)") "template", ii1, ".sco"
 665        endif
 666
 667        open (7,file=plik4)
 668
 669
 670       if (peptide.eq.1) then
 671       do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee)
 672         if (glys1.eq.0) then
 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"
 678         endif
 679
 680         if (glys1.eq.1) then
 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"
 686         endif
 687       enddo
 688       endif
 689
 690
 691       if (peptide.eq.0) then
 692       do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee)
 693 c        write(*,*) "ijl"
 694 c no Gly
 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"
 706         endif
 707
 708
 709 c chain 1 first Gly
 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"
 721         endif
 722
 723 c chain 1 last Gly
 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"
 735         endif
 736
 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"
 749         endif
 750
 751 c chain 2 first Gly
 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"
 763         endif
 764
 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"
 777         endif
 778
 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"
 791         endif
 792
 793 c all Gly
 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"
 805         endif
 806
 807
 808
 809       enddo
 810       endif
 811
 812
 813        close (7)
 814
 815 c     koniec petli po ilosci struktur
 816       enddo
 817
 818 c      do i=1,20
 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
 821 c      enddo
 822
 823 c      do i=1,20
 824 c      do j=1,20
 825 c         dd=sqrt((cmxx(i)-cmxx(j))**2+(cmyy(i)-cmyy(j))**2+
 826 c     &      (cmzz(i)-cmzz(j))**2)
 827 c         write(*,*) i,j,dd
 828 c      enddo
 829 c      enddo
 830
 831
 832   300 format(6f8.2)
 833   301 format(3f8.2)
 834   302 format(f8.2)
 835   800 format(a4,2X,i5,1X,a4,1X,a3,1X,a1,i4,4X,3f8.3)
 836   303 format(a3)
 837   304 format(a6)
 838   305 format(a5,1X,i8)
 839       end
 840
 841
 842
 843
 844       SUBROUTINE SRAND(ISEED)
 845 C
 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.
 850 C
 851       COMMON /SEED/JSEED,IFRST
 852 C
 853       JSEED = ISEED
 854       IFRST = 0
 855 C
 856       RETURN
 857       END
 858       REAL FUNCTION RAND()
 859 C
 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:
 863 C
 864 C     Park, Steven K. and Miller, Keith W., "Random Number Generators:
 865 C     Good Ones are Hard to Find", Communications of the ACM,
 866 C     October, 1988.
 867 C
 868       PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773,
 869      +           MOMDMP=2836)
 870 C
 871       COMMON  /SEED/JSEED,IFRST
 872       INTEGER HVLUE, LVLUE, TESTV, NEXTN
 873       SAVE    NEXTN
 874 C
 875       IF (IFRST .EQ. 0) THEN
 876         NEXTN = JSEED
 877         IFRST = 1
 878       ENDIF
 879 C
 880       HVLUE = NEXTN / MOBYMP
 881       LVLUE = MOD(NEXTN, MOBYMP)
 882       TESTV = MPLIER*LVLUE - MOMDMP*HVLUE
 883       IF (TESTV .GT. 0) THEN
 884         NEXTN = TESTV
 885       ELSE
 886         NEXTN = TESTV + MODLUS
 887       ENDIF
 888       RAND = REAL(NEXTN)/REAL(MODLUS)
 889 C
 890       RETURN
 891       END
 892       BLOCKDATA RANDBD
 893       COMMON /SEED/JSEED,IFRST
 894 C
 895       DATA JSEED,IFRST/123456789,0/
 896 C
 897       END
 898
 899       Subroutine ranorND(b,nd,ilosc_atomow2,ii1)
 900       double precision b
 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)
 908       integer seed
 909 c      write(*,*) "In Subroutine ranorND(b,nd), nd=", nd
 910 c      write(*,*) "ilosc_atomow2=", ilosc_atomow2
 911 c      do j=1,ilosc_atomow2
 912         ransq=2
 913         seed = time()*ii1
 914         call srand(seed)
 915         do while (ransq.ge.1)
 916 c       write(*,*) ransq
 917            rsq=0.0
 918            do i=1,nd-1
 919               ran(i)=1.0-2.0*rand(0)
 920 c           write(*,*) i, ran(i)
 921               rsq=ran(i)*ran(i)+rsq          !  (3)
 922 c           write(*,*) i, rsq
 923            enddo
 924            ransq=rsq
 925 c           write(*,*) j, ransq
 926         enddo
 927         ranh=2.0*sqrt(1.0-ransq)
 928 c       write(*,*) j, ranh, ransq
 929         do i=1,nd-1
 930           b(i)=ran(i)*ranh
 931 c          write(*,*) ran(i), ranh
 932         enddo
 933         b(nd)=(1.0-2.0*ransq)
 934 c        write(*,*) b(1,j), b(2,j), b(3,j), b(4,j)
 935 c      enddo
 936       return
 937       end
 938
 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)
 943       integer i
 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))
 954
 955 c      write(*,*) r11, r22, r33, r12, r21, r13, r31, r23, r32
 956
 957       do i=1,ilosc_atomow2
 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)
 963       enddo
 964
 965       return
 966       end
 967