From: Cezary Czaplewski Date: Wed, 26 Feb 2020 23:05:20 +0000 (+0100) Subject: unres-dock generator cmake X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=8426fef1ec5dd059f15333c150da903b74e7af69 unres-dock generator cmake --- diff --git a/CMakeLists.txt b/CMakeLists.txt index 821c630..e80732e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -212,5 +212,6 @@ else() add_subdirectory(source/cluster/unres/src) add_subdirectory(source/xdrfpdb/src) add_subdirectory(source/xdrfpdb/src-M) + add_subdirectory(source/unres-dock) endif(UNRES_NA_MMCE) diff --git a/source/unres-dock/CMakeLists.txt b/source/unres-dock/CMakeLists.txt new file mode 100644 index 0000000..f7a596f --- /dev/null +++ b/source/unres-dock/CMakeLists.txt @@ -0,0 +1,25 @@ +# + +set(UNRES_DOCK + generator_v13b.F +) + + +if (Fortran_COMPILER_NAME STREQUAL "gfortran") + set(CPPFLAGS "${CPPFLAGS} -DGFORTRAN") +endif (Fortran_COMPILER_NAME STREQUAL "gfortran") + + +#========================================= +# Build the binaries +#========================================= +add_executable(UNRES_DOCK_BIN ${UNRES_DOCK} ) +set_target_properties(UNRES_DOCK_BIN PROPERTIES OUTPUT_NAME generator ) +set_property(TARGET UNRES_DOCK_BIN PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin ) + + +#========================================= +# Install Path +#========================================= +install(TARGETS UNRES_DOCK_BIN DESTINATION ${CMAKE_INSTALL_PREFIX}/unres-dock) + diff --git a/source/unres-dock/generator_v13b.F b/source/unres-dock/generator_v13b.F new file mode 100644 index 0000000..0a33cab --- /dev/null +++ b/source/unres-dock/generator_v13b.F @@ -0,0 +1,967 @@ + program generator + implicit none + integer proba, conf, cm1, cm2, accept, nd + integer krok, i, n, j, g, ii1, ilosc, ijl + integer nnn, seed, proba1, proba2, proba3, proba0x,proba0y,proba0z + integer ilosc_atomow1, ilosc_atomow2, korekta1, korekta2 + integer ile_reszt_mol1, ile_reszt_mol2, ir, peptide,numerr2(90000) + integer numera1(90000),numerr1(90000),numera2(90000) + integer glys1, glys2, glye1, glye2, ilee + real*8 rand, r, boxx, boxy, boxz, dist2acm, dist1acm, b(4), t(3) + real*8 minx1, miny1, minz1, maxx1, maxy1, maxz1 + real*8 minx2, miny2, minz2, maxx2, maxy2, maxz2 + real*8 distx1, disty1, distz1, distx2, disty2, distz2 + real*8 xcm1, ycm1, zcm1, xcm2, ycm2, zcm2 + real*8 dist1cm, dist1cmt, dist2cm, dist2cmt + real*8 cm2xr, cm2yr, cm2zr,x2r(90000),y2r(90000),z2r(90000),dist12 + real*8 angle, c, s, x2new, y2new, z2new + character*60 plik1, plik2, plik3, plik4 + character*54 linia + real*8 x1(90000),y1(90000),z1(90000),x2(90000),y2(90000),z2(90000) + real*8 dd, wymiar, cmxx(900), cmyy(900), cmzz(900),xt2(90000,3) + real*8 xf(90000,3) + character*1 chain + character*3 r1s, r2s, r1e, r2e + character*5 atoma1(90000), atoma2(90000) + character*6 nazwaa1(90000), nazwar1(90000) + character*6 nazwaa2(90000), nazwar2(90000) + integer time_a(3) + + korekta1 = 1 + korekta2 = 1 + +c write(*,*)"Sposob uruchamiania:" +c write(*,*)"nazwa_programu nazwa_pliku_inp_pdb1 +c & nazwa_pliku_inp_pdb2 nazwa_pliku_out" +c write(*,*)"" + + nnn = iargc() + if (nnn.lt.3) then + write(*,*)"program input_pdb1 input_pdb2 + & [0-1] 0 for protein, 1 for peptide" + stop + endif + + call getarg(1,plik1) + call getarg(2,plik2) + call getarg(3,plik3) + open(1,file=plik1,status="old") + open(2,file=plik2,status="old") + open(4,file='temp1') + open(5,file='temp2') + write(*,*) "All files read" + + read(plik3(1:1),'(i1)') peptide + + if (peptide.ne.0.and.peptide.ne.1) then + write(*,*) "Third parameter must be 0 or 1, + & where 0 is for a protein-protein docking, + & and 1 is for protein-peptide docking." + STOP + endif + + write(*,*) "Option", peptide, "selected + & (0 for protein, 1 for peptide)" + + +c MOLEKULA 1 + do while (.true.) + read(1,'(a54)',end=22) linia + if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then + write(4,'(a54)') linia + endif + if (linia(1:3).eq.'TER') then +c write(2,'(a54)') linia + goto 22 + endif + enddo + +22 continue + close(4) + write(*,*) "Molecule 1 read" + + +c MOLEKULA 2 + do while (.true.) + read(2,'(a54)',end=23) linia + if (linia(1:4).eq.'ATOM'.and.linia(17:17).lt.'B') then + write(5,'(a54)') linia + endif + if (linia(1:3).eq.'TER') then +c write(2,'(a54)') linia + goto 23 + endif + enddo +23 continue + close(5) + write(*,*) "Molecule 2 read" + + open(4,file='temp1') + open(5,file='temp2') + + minx1=999.99 + miny1=999.99 + minz1=999.99 + + maxx1=-999.99 + maxy1=-999.99 + maxz1=-999.99 + + + minx2=999.99 + miny2=999.99 + minz2=999.99 + + maxx2=-999.99 + maxy2=-999.99 + maxz2=-999.99 + + + +c MOLEKULA 1 + i=1 + xcm1=0 + ycm1=0 + zcm1=0 + ile_reszt_mol1=0 + + do while (.true.) + read(4,800,end=111) + & atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),chain,numerr1(i), + & x1(i),y1(i),z1(i) + + +c write(*,300) minx, maxx, miny, maxy, minz, maxz + + if (x1(i).lt.minx1) minx1=x1(i) + if (y1(i).lt.miny1) miny1=y1(i) + if (z1(i).lt.minz1) minz1=z1(i) + + if (x1(i).gt.maxx1) maxx1=x1(i) + if (y1(i).gt.maxy1) maxy1=y1(i) + if (z1(i).gt.maxz1) maxz1=z1(i) + +c write(*,300) minx, maxx, miny, maxy, minz, maxz + +c write(*,*) nazwaa1(i) + if (nazwaa1(i).eq.' CA ') then + xcm1 = xcm1 + x1(i) + ycm1 = ycm1 + y1(i) + zcm1 = zcm1 + z1(i) + ile_reszt_mol1 = ile_reszt_mol1 +1 + endif + + if (atoma1(i).eq.'TER ') goto 111 + if (atoma1(i).eq.'END ') goto 111 + if (i.eq.1) write(3,'(a6)') "REMARK" +c write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i), +c & numerr1(i),x1(i),y1(i),z1(i) + +c write(*,*) i, x(i),y(i),z(i) + + i=i+1 + enddo +111 continue + + distx1=maxx1-minx1 + disty1=maxy1-miny1 + distz1=maxz1-minz1 + ilosc_atomow1 = i-1 + xcm1 = xcm1 / (ile_reszt_mol1*1.0) + ycm1 = ycm1 / (ile_reszt_mol1*1.0) + zcm1 = zcm1 / (ile_reszt_mol1*1.0) + + do i=1,ilosc_atomow1 + x1(i)=x1(i)-xcm1 + y1(i)=y1(i)-ycm1 + z1(i)=z1(i)-zcm1 + enddo + + xcm1=0 + ycm1=0 + zcm1=0 + + do i=1,ilosc_atomow1 + if (nazwaa1(i).eq.' CA ') then + xcm1 = xcm1 + x1(i) + ycm1 = ycm1 + y1(i) + zcm1 = zcm1 + z1(i) + endif + enddo + xcm1 = xcm1 / (ile_reszt_mol1*1.0) + ycm1 = ycm1 / (ile_reszt_mol1*1.0) + zcm1 = zcm1 / (ile_reszt_mol1*1.0) + + dist1cm=999.9 + dist1acm=0.0 + do i=1,ilosc_atomow1 +c write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i), +c & numerr1(i),x1(i)-xcm1,y1(i)-ycm1,z1(i)-zcm1 + + dist1cmt = dsqrt(((xcm1-x1(i))**2)+ + & ((ycm1-y1(i))**2)+ + & ((zcm1-z1(i))**2)) + +c write(*,*) dist1cmt, dist1cm, i, cm1 + if (dist1cm.gt.dist1cmt) then + dist1cm=dist1cmt + cm1=i + endif + + enddo + + do i=1,ilosc_atomow1 + do j=i,ilosc_atomow1 + dist1cmt = dsqrt(((x1(j)-x1(i))**2)+ + & ((y1(j)-y1(i))**2)+ + & ((z1(j)-z1(i))**2)) + if (dist1acm.lt.dist1cmt) then + dist1acm=dist1cmt + endif + enddo + enddo + + +c write(3,303) "TER" +c write(3,304) "ENDMDL" + + write(*,*) "Molecule 1 calculated" + +c MOLEKULA 2 + i=1 + xcm2=0 + ycm2=0 + zcm2=0 + ile_reszt_mol2=0 + + do while (.true.) + read(5,800,end=112) + & atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),chain,numerr2(i), + & x2(i),y2(i),z2(i) + +c write(*,300) minx, maxx, miny, maxy, minz, maxz + + if (x2(i).lt.minx2) minx2=x2(i) + if (y2(i).lt.miny2) miny2=y2(i) + if (z2(i).lt.minz2) minz2=z2(i) + + if (x2(i).gt.maxx2) maxx2=x2(i) + if (y2(i).gt.maxy2) maxy2=y2(i) + if (z2(i).gt.maxz2) maxz2=z2(i) + +c write(*,300) minx, maxx, miny, maxy, minz, maxz + +c write(*,*) nazwaa1(i) + if (nazwaa2(i).eq.' CA ') then + xcm2 = xcm2 + x2(i) + ycm2 = ycm2 + y2(i) + zcm2 = zcm2 + z2(i) + ile_reszt_mol2 = ile_reszt_mol2 + 1 + endif + + if (atoma2(i).eq.'TER ') goto 111 + if (atoma2(i).eq.'END ') goto 111 +c if (i.eq.1) write(3,'(a6)') "REMARK" +c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i), +c & numerr2(i),x2(i),y2(i),z2(i) +c write(*,*) i, x(i),y(i),z(i) + + i=i+1 + enddo +112 continue + + distx2=maxx2-minx2 + disty2=maxy2-miny2 + distz2=maxz2-minz2 + ilosc_atomow2 = i-1 + + xcm2 = xcm2 / (ile_reszt_mol2*1.0) + ycm2 = ycm2 / (ile_reszt_mol2*1.0) + zcm2 = zcm2 / (ile_reszt_mol2*1.0) + + do i=1,ilosc_atomow2 + x2(i)=x2(i)-xcm2 + y2(i)=y2(i)-ycm2 + z2(i)=z2(i)-zcm2 +c write(*,*) x2(i),y2(i),z2(i) + xt2(i,1)=x2(i) + xt2(i,2)=y2(i) + xt2(i,3)=z2(i) +c write(*,*) xt2(i,1), xt2(i,2), xt2(i,3) + enddo + + xcm2=0 + ycm2=0 + zcm2=0 + do i=1,ilosc_atomow2 + if (nazwaa2(i).eq.' CA ') then + xcm2 = xcm2 + x2(i) + ycm2 = ycm2 + y2(i) + zcm2 = zcm2 + z2(i) + endif + enddo + xcm2 = xcm2 / (ile_reszt_mol2*1.0) + ycm2 = ycm2 / (ile_reszt_mol2*1.0) + zcm2 = zcm2 / (ile_reszt_mol2*1.0) + + + dist2cm=999.9 + dist2acm=0.0 + do i=1,ilosc_atomow2 +c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i), +c & numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2 + + dist2cmt = dsqrt(((xcm2-x2(i))**2)+ + & ((ycm2-y2(i))**2)+ + & ((zcm2-z2(i))**2)) + +c write(*,*) dist2cmt, dist2cm, i, cm2 + if (dist2cm.gt.dist2cmt) then + dist2cm=dist2cmt + cm2=i + endif + + enddo + + + do i=1,ilosc_atomow2 + do j=i,ilosc_atomow2 + dist2cmt = dsqrt(((x2(j)-x2(i))**2)+ + & ((y2(j)-y2(i))**2)+ + & ((z2(j)-z2(i))**2)) + if (dist2acm.lt.dist2cmt) then + dist2acm=dist2cmt + endif + enddo + enddo + + +c write(3,303) "TER" + + glys1=0 + glys2=0 + glye1=0 + glye2=0 + r1s=nazwar1(1) + r2s=nazwar2(1) + r1e=nazwar1(ile_reszt_mol1) + r2e=nazwar2(ile_reszt_mol2) + + if (r1s.eq."GLY") then + glys1=1 + write(*,*) "Chain 1 is starting from Gly ", r1s, glys1 + endif + + if (r2s.eq."GLY") then + glys2=1 + write(*,*) "Chain 2 is starting from Gly ", r2s, glys2 + endif + + if (r1e.eq."GLY") then + glye1=1 + write(*,*) "Chain 1 is ending at Gly ", r1e, glye1 + endif + + + if (r2e.eq."GLY") then + glye2=1 + write(*,*) "Chain 2 is ending from Gly ", r2e, glye2 + endif + + +c checking for Gly + ilee=4 + if(glys1.eq.1) ilee=ilee-1 + if(glys2.eq.1) ilee=ilee-1 + if(glye1.eq.1) ilee=ilee-1 + if(glye2.eq.1) ilee=ilee-1 + write(*,*) "How many dummy atoms:", ilee + write(*,*) "" + + + write(*,*) "Molecule 2 calculated" + write(*,*) "Number of atoms:", ilosc_atomow1, ilosc_atomow2 + + write(*,*) "Min and max coordinates in the systems" + write(*,300) minx1, maxx1, miny1, maxy1, minz1, maxz1 + write(*,300) minx2, maxx2, miny2, maxy2, minz2, maxz2 + + write(*,*) "Distances in axes x, y, z" + write(*,301) distx1, disty1, distz1 + write(*,301) distx2, disty2, distz2 + + boxx = (distx1 + distx2)*1.5 + 20.0 + boxy = (disty1 + disty2)*1.5 + 20.0 + boxz = (distz1 + distz2)*1.5 + 20.0 + +c write(*,*) "Boxsize [x, y, z]" +c write(*,301) boxx, boxy, boxz + + write(*,*) "Maximum distances:" + write(*,*) dist1acm, dist2acm + + write(*,*) "Boxsize (suggested)" + write(*,*) (dist1acm+dist2acm)*1.2+20 + open(61,file='boxsize.txt') + write(61,302) ((dist1acm+dist2acm)*1.2+20) + close (61) + + write(*,*) "Center of masses:" + write(*,301) xcm1, ycm1, zcm1 + write(*,301) xcm2, ycm2, zcm2 + + write(*,*) "Center of masses atoms:" + write(*,800) atoma1(cm1),numera1(cm1),nazwaa1(cm1),nazwar1(cm1), + & " ",numerr1(cm1),x1(cm1),y1(cm1),z1(cm1) + write(*,800) atoma2(cm2),numera2(cm2),nazwaa2(cm2),nazwar2(cm2), + & " ",numerr2(cm2),x2(cm2),y2(cm2),z2(cm2) + + + write(*,*) "Number of amino-acid residues in two proteins:" + write(*,*) ile_reszt_mol1, ile_reszt_mol2 + + +c --------------------------------------------------------------------- +c KONIEC WCZYTYWANIA +c --------------------------------------------------------------------- + + korekta1 = numerr1(1) + write(*,*) "Original number of the first residue in + & the first protein:" + write(*,*) numerr1(1) + + write(*,*) "Correction of the first residue in the first protein:" + write(*,*) korekta1 + + +c wymiar=(distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2) +c wymiar=wymiar/12.0 +c wymiar=10.0 + +c write(*,*) "Minimum required distance between centers of mass:" +c write(*,*) wymiar + + boxx = (7+distx1 + dist2acm)/2.0 + boxy = (7+disty1 + dist2acm)/2.0 + boxz = (7+distz1 + dist2acm)/2.0 + write(*,*) "Rotation space (+/-):", boxx, boxy, boxz + + + do ii1=1,24 + + proba = 1 + proba0x = 1 + proba0y = 1 + proba0z = 1 + proba1 = 1 + proba2 = 1 + proba3 = 1 + conf = 1 + + + if (ii1.lt.10) then + write (plik3, "(A6,I1,a4)") "model0", ii1, ".pdb" + else + write (plik3, "(A5,I2,a4)") "model", ii1, ".pdb" + endif + + write(*,*) "Conformation will be saved to file:", plik3 + +c plik3a = plik3//ilosc + open(3,file=plik3) + + + do i=1,ilosc_atomow1 + write(3,800) atoma1(i),numera1(i),nazwaa1(i),nazwar1(i),"A", + & (numerr1(i)-korekta1+1), + & x1(i),y1(i),z1(i) + enddo + write(3,303) "TER" +c write(3,304) "ENDMDL" + + if (peptide.eq.1) then + wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/4 + endif + + if (peptide.eq.0) then + wymiar=((distx1 + distx2)+(disty1 + disty2)+(distz1 + distz2))/7 + endif + + if (ii1.eq.1) then + write(*,*) "Initial minimum distance between centers of mass:", + & wymiar + endif + +c ROTATIONS OF THE SECOND MOLECULE + nd=4 + t(1)=0.0 + t(2)=0.0 + t(3)=0.0 + call ranorND(b,nd,ilosc_atomow2,ii1) +c do i=1,4 +c write(*,*) b(i) +c enddo + call rot_trans(xt2,t,b,xf,ilosc_atomow2) + do i=1,ilosc_atomow2 + x2(i)=xf(i,1) + y2(i)=xf(i,2) + z2(i)=xf(i,3) + enddo +c END OF ROTATIONS + + + do while (.true.) + if (conf.gt.1) goto 133 + +#ifdef GFORTRAN + seed = time()*proba +#else + call itime(time_a) + seed = (60*time_a(2)+time_a(3))*proba +#endif + call srand(seed) + r = rand(0) +c write(*,*)rand(0),rand(0),rand(0) +c write(*,*)rand(0),rand(0),rand(0) + + cm2xr = (rand(0)*2*boxx)-(boxx) + cm2yr = (rand(0)*2*boxy)-(boxy) + cm2zr = (rand(0)*2*boxz)-(boxz) +c write(*,*) cm2xr, cm2yr, cm2zr +c write(*,300)rand(0),rand(0),rand(0) + +c write(*,*) ABS(cm2xr+xcm1), (dist1acm+dist2acm+3) + +c if (ABS(cm2xr+xcm1).gt.((dist1acm+dist2acm+5)*0.8)) then +c proba0x=proba0x+1 +c goto 200 +c endif + +c if (ABS(cm2yr+ycm1).gt.((dist1acm+dist2acm+5)*0.8)) then +c proba0y=proba0y+1 +c goto 200 +c endif + +c if (ABS(cm2zr+zcm1).gt.((dist1acm+dist2acm+5)*0.8)) then +c proba0z=proba0z+1 +c goto 200 +c endif + + + do i=1,ilosc_atomow2 + x2r(i)=x2(i)-cm2xr + y2r(i)=y2(i)-cm2yr + z2r(i)=z2(i)-cm2zr + enddo + + +c SPRAWDZANIE CZY NIE JEST ZBYT PODOBNE DO JUZ WYLOSOWANEGO + cmxx(ii1)=cm2xr + cmyy(ii1)=cm2yr + cmzz(ii1)=cm2zr + + if (ii1.gt.1) then + + do ir=1,ii1-1 + dd=sqrt((cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+ + & (cm2zr-cmzz(ir))**2) + +c write(*,*) "Center of mass is", dd, "A away from", ir +c write(*,*) "Minimum distance required:", wymiar, "A" + +c write(*,*) ii1, ir, proba3, dd, wymiar + + if (dd.lt.wymiar) then + proba3=proba3+1 + if (proba3.gt.1.and.peptide.eq.1) wymiar=wymiar*0.99992 + if (proba3.gt.1.and.peptide.eq.0) wymiar=wymiar*0.99992 +c write(*,*) proba3, wymiar + goto 200 + endif + + enddo + endif + + + + +c SPRAWDZANIE ODLEGLOSCI + accept=0 + do i=1,ilosc_atomow2 +c write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i), +c & numerr2(i),x2(i)-xcm2,y2(i)-ycm2,z2(i)-zcm2 + + do j=1,ilosc_atomow1 + dist12 = sqrt(((x2r(i)-x1(j))**2)+ + & ((y2r(i)-y1(j))**2)+ + & ((z2r(i)-z1(j))**2)) +c write(*,*) dist12, cm2xr, cm2yr, cm2zr + if (dist12.lt.2.5) then + proba1=proba1+1 + goto 200 + endif + + if (dist12.le.3.5) accept=1 +c if (dist12.le.3.5) accept=1 +c if (dist12.gt.3.5) proba2=proba2+1 + enddo + enddo + +c write(*,*) "No overlaps", proba + if (accept.eq.0) goto 200 +c write(*,*) "No overlaps, contact", conf, proba +c write(3,305) "MODEL", conf +c conf = conf+1 + + + + write(*,*) "Conformation", ii1, "Is accepted after", + & proba0x+proba0y+proba0z, + & proba1, proba2, proba3, proba, "attempts. ", + & "Dist between CoM of previous=", wymiar + +c Dodatkowe info + write(*,*)"Extr",ABS(cm2xr+xcm1), ABS(cm2yr+ycm1), ABS(cm2zr+zcm1) + + + conf=conf+1 +c write (*,*) + + + + + + korekta2 = numerr2(1)-1 +c write(*,*) "Original number of the first residue in +c & the second protein:" +c write(*,*) numerr2(1) + +c write(*,*)"Correction of the first residue in the second protein:" +c write(*,*) korekta2 + do i=1,ilosc_atomow2 + write(3,800) atoma2(i),numera2(i),nazwaa2(i),nazwar2(i),"B", + & (numerr2(i)+ile_reszt_mol1-korekta2), + & x2r(i),y2r(i),z2r(i) + enddo + write(3,303) "TER" +c write(3,304) "ENDMDL" +c write(3,305) "MODEL", conf + write(3,303) "END" + + 200 continue + proba = proba + 1 + enddo + 133 continue + + close (3) + +c call sleep(1) + + + if (ii1.lt.10) then + write (plik4, "(A9,I1,a4)") "template0", ii1, ".sco" + else + write (plik4, "(A8,I2,a4)") "template", ii1, ".sco" + endif + + open (7,file=plik4) + + + if (peptide.eq.1) then + do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee) + if (glys1.eq.0) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + + if (glys1.eq.1) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+1)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + enddo + endif + + + if (peptide.eq.0) then + do ijl=1,(ile_reszt_mol1+ile_reszt_mol2+ilee) +c write(*,*) "ijl" +c no Gly + if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.0) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1+3).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+3)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" + if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+4)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + + +c chain 1 first Gly + if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.0) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1+2).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" + if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + +c chain 1 last Gly + if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.0) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1+2).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" + if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + +c chain 1 first Gly and chain 1 last Gly + if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.0) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1+1).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" + if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + +c chain 2 first Gly + if (glys1.eq.0.and.glye1.eq.0.and.glys2.eq.1) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1+2).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" + if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+3)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + +c chain 1 first Gly chain 2 first Gly + if (glys1.eq.1.and.glye1.eq.0.and.glys2.eq.1) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.ge.(ile_reszt_mol1+1).and.ijl.le.(ile_reszt_mol1+1)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1+1).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" + if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + +c chain 1 last Gly chain 2 first Gly + if (glys1.eq.0.and.glye1.eq.1.and.glys2.eq.1) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.le.ile_reszt_mol1+1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" +c if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3)) +c & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1+1).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1+1)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" + if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+2)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + +c all Gly + if (glys1.eq.1.and.glye1.eq.1.and.glys2.eq.1) then + if (ijl.eq.1) write(7,*) ijl, "1.0 1.0 1.0 ", "1" + if (ijl.le.ile_reszt_mol1.and.ijl.gt.1) + & write(7,*) ijl, "1.0 1.0 1.0 ", "1" +c if (ijl.ge.(ile_reszt_mol1+2).and.ijl.le.(ile_reszt_mol1+3)) +c & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + if (ijl.gt.(ile_reszt_mol1).and. + & ijl.le.(ile_reszt_mol2+ile_reszt_mol1)) + & write(7,*) ijl, "1.0 1.0 1.0 ", "2" +c if (ijl.ge.(ile_reszt_mol2+ile_reszt_mol1+1)) +c & write(7,*) ijl, "1.0 1.0 1.0 ", "0" + endif + + + + enddo + endif + + + close (7) + +c koniec petli po ilosci struktur + enddo + +c do i=1,20 +c write(*,*) i, cmxx(i), cmyy(i), cmzz(i) +c dd=(cm2xr-cmxx(ir))**2+(cm2yr-cmyy(ir))**2+(cm2zr-cmyy(ir))**2 +c enddo + +c do i=1,20 +c do j=1,20 +c dd=sqrt((cmxx(i)-cmxx(j))**2+(cmyy(i)-cmyy(j))**2+ +c & (cmzz(i)-cmzz(j))**2) +c write(*,*) i,j,dd +c enddo +c enddo + + + 300 format(6f8.2) + 301 format(3f8.2) + 302 format(f8.2) + 800 format(a4,2X,i5,1X,a4,1X,a3,1X,a1,i4,4X,3f8.3) + 303 format(a3) + 304 format(a6) + 305 format(a5,1X,i8) + end + + + + + SUBROUTINE SRAND(ISEED) +C +C This subroutine sets the integer seed to be used with the +C companion RAND function to the value of ISEED. A flag is +C set to indicate that the sequence of pseudo-random numbers +C for the specified seed should start from the beginning. +C + COMMON /SEED/JSEED,IFRST +C + JSEED = ISEED + IFRST = 0 +C + RETURN + END + REAL FUNCTION RAND() +C +C This function returns a pseudo-random number for each invocation. +C It is a FORTRAN 77 adaptation of the "Integer Version 2" minimal +C standard number generator whose Pascal code appears in the article: +C +C Park, Steven K. and Miller, Keith W., "Random Number Generators: +C Good Ones are Hard to Find", Communications of the ACM, +C October, 1988. +C + PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773, + + MOMDMP=2836) +C + COMMON /SEED/JSEED,IFRST + INTEGER HVLUE, LVLUE, TESTV, NEXTN + SAVE NEXTN +C + IF (IFRST .EQ. 0) THEN + NEXTN = JSEED + IFRST = 1 + ENDIF +C + HVLUE = NEXTN / MOBYMP + LVLUE = MOD(NEXTN, MOBYMP) + TESTV = MPLIER*LVLUE - MOMDMP*HVLUE + IF (TESTV .GT. 0) THEN + NEXTN = TESTV + ELSE + NEXTN = TESTV + MODLUS + ENDIF + RAND = REAL(NEXTN)/REAL(MODLUS) +C + RETURN + END + BLOCKDATA RANDBD + COMMON /SEED/JSEED,IFRST +C + DATA JSEED,IFRST/123456789,0/ +C + END + + Subroutine ranorND(b,nd,ilosc_atomow2,ii1) + double precision b +C generates a random vector on a unit sphere onC an ND-dimensional space. +C algorithm 40 page 410. +C the algorithm assumes that the formula for bz is valid, i.e +C bz=(1.0-2.0*ransq), +C and also that, ransq can be computed as as the sum of the independent random +C numbers to the second power (expression 3). + dimension ran(nd-1), b(nd) + integer seed +c write(*,*) "In Subroutine ranorND(b,nd), nd=", nd +c write(*,*) "ilosc_atomow2=", ilosc_atomow2 +c do j=1,ilosc_atomow2 + ransq=2 + seed = time()*ii1 + call srand(seed) + do while (ransq.ge.1) +c write(*,*) ransq + rsq=0.0 + do i=1,nd-1 + ran(i)=1.0-2.0*rand(0) +c write(*,*) i, ran(i) + rsq=ran(i)*ran(i)+rsq ! (3) +c write(*,*) i, rsq + enddo + ransq=rsq +c write(*,*) j, ransq + enddo + ranh=2.0*sqrt(1.0-ransq) +c write(*,*) j, ranh, ransq + do i=1,nd-1 + b(i)=ran(i)*ranh +c write(*,*) ran(i), ranh + enddo + b(nd)=(1.0-2.0*ransq) +c write(*,*) b(1,j), b(2,j), b(3,j), b(4,j) +c enddo + return + end + + subroutine rot_trans(xt2,t,q,xf,ilosc_atomow2) +C x coordinates, t translation vector, q quaternion + double precision xt2(90000,3),xf(90000,3) + double precision t(3),q(4) + integer i + write(*,*) "Rotation vector:", q(1), q(2), q(3), q(4) + r11= q(1)*q(1)+q(2)*q(2)-q(3)*q(3)-q(4)*q(4) + r22= q(1)*q(1)-q(2)*q(2)+q(3)*q(3)-q(4)*q(4) + r33= q(1)*q(1)-q(2)*q(2)-q(3)*q(3)+q(4)*q(4) + r12= 2*(q(2)*q(3)-q(1)*q(4)) + r21= 2*(q(2)*q(3)+q(1)*q(4)) + r13= 2*(q(2)*q(4)+q(1)*q(3)) + r31= 2*(q(2)*q(4)-q(1)*q(3)) + r23= 2*(q(3)*q(4)-q(1)*q(2)) + r32= 2*(q(3)*q(4)+q(1)*q(2)) + +c write(*,*) r11, r22, r33, r12, r21, r13, r31, r23, r32 + + do i=1,ilosc_atomow2 + xf(i,1)= r11*xt2(i,1)+r12*xt2(i,2)+r13*xt2(i,3) + t(1) + xf(i,2)= r21*xt2(i,1)+r22*xt2(i,2)+r23*xt2(i,3) + t(2) + xf(i,3)= r31*xt2(i,1)+r32*xt2(i,2)+r33*xt2(i,3) + t(3) +c write(*,*) i, xt2(i,1), xt2(i,2), xt2(i,3) +c write(*,*) i, xf(i,1), xf(i,2), xf(i,3) + enddo + + return + end +