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