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