! Projet de simulation numérique ! student : Christian Schlatter, SB, section physique, EPFL ! professeur : John Maddocks, SB, section mathématique, EPFL ! This release meets the version 2.2 of the PDB format description guide ! which is available on ! http://www.rcsb.org/pdb/docs/format/pdbguide2.2/guide2.2_frame.html ! Future modifications of the field lengths (as proposed for the atom entry) ! will require a modification of this code. module parameters ! User adjustable parameter module integer, parameter :: memory = 12 ! Memory limit of smallest radii integer, parameter :: limit = 200 ! Upper limit for distribution vector logical, parameter :: talkmode = .false. ! Comments on screen logical, parameter :: secondrun = .true. ! Statistics run to classify distribution of radii character (LEN=25), parameter :: releaseinfo = 'Release 1.61 (2002/07/07)' end module parameters ! ********************************************************************** module var ! Data type definitions type carbons integer :: residue real, dimension(3) :: coordinates character (LEN=1) :: chainident end type carbons end module var ! ********************************************************************** module mathe ! Mathematical routines contains function crossproduct (a, b) ! Calculates the cross product of two vectors implicit none real, dimension(3), intent(in) :: a, b real, dimension(3) :: crossproduct crossproduct(1) = a(2)*b(3) - a(3)*b(2) crossproduct(2) = a(3)*b(1) - a(1)*b(3) crossproduct(3) = a(1)*b(2) - a(2)*b(1) end function crossproduct function determinant (A) ! Calculates the determinant of a matrix using ! the partial pivoting Gaussian elimination scheme. implicit none real, pointer, dimension(:,:) :: A real, allocatable, dimension(:,:) :: B real, allocatable, dimension(:) :: C integer, allocatable, dimension(:) :: indexer integer :: i, j, k, n, tempindex, msgn real :: scaler, pivot, pivoti, pivotj real :: determinant n = size(A(:,1)) allocate (B(n,n)) B = A ! index records the pivoting order allocate (indexer(n)) ! A is the original matrix in the input. B is A transformed ! into the matrix plus the pivoting element ratios below ! the diagonal in the output. do i = 1, n indexer(i) = i end do ! Find the rescaling factors, one from each row allocate (C(n)) do i = 1, n scaler = 0.0 do j = 1, n scaler = amax1(scaler,abs(B(i,j))) end do C(i) = scaler end do ! Search the pivoting (largest) element from each column do j = 1, n-1 pivot = 0.0 do i = j, n pivoti = abs(B(indexer(i),j)) / C(indexer(i)) if (pivoti.GT.pivot) then pivot = pivoti k = i endif end do ! Interchange the rows via indexer(n) to record pivoting order tempindex = indexer(j) indexer(j) = indexer(k) indexer(k) = tempindex do i = j+1, n pivotj = B(indexer(i),j)/B(indexer(j),j) ! Record pivoting ratios below the diagonal B(indexer(i),j) = pivotj ! Modify other elements accordingly do k = j+1, n B(indexer(i),k) = B(indexer(i),k)-pivotj*B(indexer(j),k) end do end do end do determinant = 1.0 do i = 1, n determinant = determinant*B(indexer(i),i) end do msgn = 1 do i = 1, n do while (i.NE.indexer(i)) msgn = -msgn j = indexer(i) indexer(i) = indexer(j) indexer(j) = j end do end do determinant = msgn*determinant deallocate(indexer) deallocate(B) deallocate(C) end function determinant end module mathe ! ********************************************************************** module distance implicit none contains subroutine circles (alphacarbons, f, g, truelocal, globalradius, & globallocation, smallestradius, smallestlocation, & distribution, thickness) ! fits circles through the macromolecule use var use mathe use parameters implicit none type(carbons), pointer, dimension(:) :: alphacarbons real :: r2, A, B, C, thickness real, dimension (3) :: D, first, second, third real, dimension (3) :: crossed, radii real, pointer, dimension(:) :: f, g, localradius, nonadjacentradius real, pointer, dimension(:) :: globalradius, nonlocalradius real, pointer, dimension(:) :: partlyadjacentradius, truelocal real, pointer, dimension(:,:) :: smallestradius integer, pointer, dimension(:) :: firstlocal integer, pointer, dimension(:,:) :: nonadjacentlocation integer, pointer, dimension(:,:) :: partlyadjacentlocation integer, pointer, dimension(:,:) :: globallocation integer, pointer, dimension(:,:,:) :: smallestlocation integer, dimension(3) :: combination integer :: numbercarbon, i, j, k, l, m, q, r integer :: mem, index, pos, location integer, dimension(1) :: selector ! due to compatibility issues with Fortran 90 integer(4), dimension(limit) :: distribution numbercarbon = size (alphacarbons) allocate (localradius(numbercarbon)) allocate (truelocal(numbercarbon)) allocate (nonadjacentradius(numbercarbon)) allocate (firstlocal(numbercarbon)) allocate (nonadjacentlocation(numbercarbon,3)) allocate (partlyadjacentradius(numbercarbon)) allocate (partlyadjacentlocation(numbercarbon,3)) if (memory .gt. numbercarbon) then ! I can't remember more than present atoms ! mem = numbercarbon else mem = memory end if localradius = 0.0 truelocal = 0.0 nonadjacentradius = 0.0 partlyadjacentradius = 0.0 smallestradius = 0.0 index = 0 if(thickness .ne. 0.0) then ! Statistical classification of the radii possible now write(*,'(" Second run. Statistics running...")') write(*, *) write(*, *) distribution = 0.0 end if ! Checking all possible combination of 3 points do i = 1 , numbercarbon - 2 do j = i + 1 , numbercarbon - 1 do k = j + 1 , numbercarbon ! colinearity check -> ((b-a)\times(c-a))^2 = 0 first = alphacarbons(i)%coordinates second = alphacarbons(j)%coordinates third = alphacarbons(k)%coordinates crossed = crossproduct((second-first), (third-first)) ! radii calculation if (dot_product(crossed, crossed) == 0) then r2 = 0.0 else ! Radius formula when three points are given D(1) = (sum((alphacarbons(i)%coordinates - & alphacarbons(j)%coordinates)**2)) D(2) = (sum((alphacarbons(i)%coordinates - & alphacarbons(k)%coordinates)**2)) D(3) = (sum((alphacarbons(j)%coordinates - & alphacarbons(k)%coordinates)**2)) A = sum(D**2) B = 2*((D(1)*D(2))+(D(1)*D(3))+(D(2)*D(3))) if((A-B) == 0) then C = -1 ! The points are colinear else C = (- product(D) / (A - B)) end if if (C < 0) then ! This may also occur if the points are almost colinear, ! so that the radius of the circle is quite large write(*,'(" Problem with alpha carbons : ", i10, ", ", i10, ", ", i10)') i, j, k write(*,'(" Colinearity check dot product equals: ", f10.3)') dot_product(crossed, crossed) write(*,*) r2 = 0.0 else r2 = sqrt (C) if(thickness .ne. 0.0) then ! This is the second run : ! Classification of the found radius pos = floor(((r2/thickness) - 1) * 100) + 1 if (pos .gt. limit) then ! Found radius is bigger than the given limit pos = limit end if distribution(pos) = distribution(pos) + 1 end if combination(1) = i ! This is due to Fortran 90 compatibility combination(2) = j combination(3) = k do m = 1, 3 ! Completing for all three points l = combination(m) do q = 1, mem if (r2 .lt. smallestradius(l,q)) then ! A smaller radius has been found do r = mem, q+1, -1 ! Shifting the radius vector for insertion of the ! new radius smallestradius(l,r) = smallestradius(l,r-1) smallestlocation(l,r,1) = smallestlocation(l,r-1,1) smallestlocation(l,r,2) = smallestlocation(l,r-1,2) smallestlocation(l,r,3) = smallestlocation(l,r-1,3) end do ! Saving the new radius smallestradius(l,q) = r2 smallestlocation(l,q,1) = i smallestlocation(l,q,2) = j smallestlocation(l,q,3) = k exit end if if (smallestradius(l,q) == 0.0) then ! This is the first radius for the given point smallestradius(l,q) = r2 smallestlocation(l,q,1) = i smallestlocation(l,q,2) = j smallestlocation(l,q,3) = k exit end if end do end do end if end if if (((k-j) + (j-i)) == 2) then ! The found radius is local combination(1) = i combination(2) = j combination(3) = k if ((truelocal(j) == 0) .or. (truelocal(j) > r2)) then ! Only for the middle point j (i < j < k) this is a true ! local radius truelocal(j) = r2 end if do m = 1, 3 l = combination(m) ! Completion for all three points if ((localradius(l) == 0) .or. (localradius(l) > r2)) then localradius(l) = r2 firstlocal(l) = i ! Indicates the first point end if end do elseif (((k-j) .gt. 1) .and. ((j-i) .gt. 1)) then ! None of the atoms are neighbors combination(1) = i combination(2) = j combination(3) = k do m = 1, 3 l = combination(m) if ((nonadjacentradius(l) == 0) .or. (nonadjacentradius(l) > r2)) then nonadjacentradius(l) = r2 nonadjacentlocation(l,1) = i nonadjacentlocation(l,2) = j nonadjacentlocation(l,3) = k end if end do else ! Two of the three atoms are direct neighbors combination(1) = i combination(2) = j combination(3) = k do m = 1, 3 l = combination(m) if ((((k-j) .ne. 1) .and. ((j-i) .ne. 1)) .or. (((k-j) == 1) .and. ((j-i) == 1))) then write (*, '("Error !!!")') end if if ((partlyadjacentradius(l) == 0) .or. (partlyadjacentradius(l) > r2)) then partlyadjacentradius(l) = r2 partlyadjacentlocation(l,1) = i partlyadjacentlocation(l,2) = j partlyadjacentlocation(l,3) = k end if end do end if index = index + 1 end do end do end do if (talkmode) then write (*, '(" ", i10, " circles have been calculated")') index write (*, *) end if ! Finalizing allocate (globalradius(numbercarbon)) allocate (globallocation(numbercarbon,3)) allocate (nonlocalradius(numbercarbon)) do i = 1, numbercarbon ! Which class of radii provides the smallest radius? radii(1) = localradius(i) radii(2) = nonadjacentradius(i) radii(3) = partlyadjacentradius(i) selector = minloc(radii, mask = radii .ne. 0.0) location = selector(1) ! Global radius lists the smallest circle for each point globalradius(i) = radii(location) select case (location) ! Indexes of the points realizing the smallest circle case (1) globallocation(i,1) = firstlocal(i) globallocation(i,2) = firstlocal(i) + 1 globallocation(i,3) = firstlocal(i) + 2 case (2) globallocation(i,:) = nonadjacentlocation(i,:) case (3) globallocation(i,:) = partlyadjacentlocation(i,:) end select ! Doing the same but by ignoring the local radii radii(1) = 0.0 radii(2) = nonadjacentradius(i) radii(3) = partlyadjacentradius(i) selector = minloc(radii, mask = radii .ne. 0.0) location = selector(1) nonlocalradius(i) = radii(location) end do allocate (f(numbercarbon)) allocate (g(numbercarbon)) f = 0.0 g = 0.0 do i = 2, numbercarbon - 1 ! Calculation of the f and g factor f(i) = nonlocalradius(i)/truelocal(i) g(i) = globalradius(i)/truelocal(i) end do deallocate (nonadjacentradius, nonadjacentlocation) deallocate (partlyadjacentradius, partlyadjacentlocation) deallocate (nonlocalradius, firstlocal, localradius) end subroutine circles ! ********************************************************************** subroutine spheres (alphacarbons, f, g, localradius, globalradius, & globallocation, smallestradius, smallestlocation, & distribution, thickness) ! fits spheres through the macromolecule use var use mathe use parameters implicit none type(carbons), pointer, dimension(:) :: alphacarbons real :: A, B, maxradius, rquadrat real :: r3, thickness, detn, detz real :: prevtime, time real, dimension(4,4) :: D, E, H real, dimension(3) :: C, first, second, third, fourth real, dimension(2) :: radii real, dimension(4) :: minima real, pointer, dimension(:) :: f, g, localradius, nonlocalradius real, pointer, dimension(:) :: globalradius real, pointer, dimension(:,:) :: smallestradius real, pointer, dimension(:,:) :: coplanarity, N, Z integer :: numbercarbon, i, j, k, l, m, o integer :: q, r, s, index, error integer :: location, pos, mem integer, pointer, dimension(:) :: firstlocal integer, pointer, dimension(:,:) :: nonlocallocation, globallocation integer, pointer, dimension(:,:,:) :: smallestlocation integer, dimension(4) :: combination integer, dimension(1) :: selector integer(4), dimension(limit) :: distribution numbercarbon = size (alphacarbons) allocate (localradius(numbercarbon)) allocate (firstlocal(numbercarbon)) allocate (nonlocalradius(numbercarbon)) allocate (nonlocallocation(numbercarbon,4)) if (memory .gt. numbercarbon) then ! I can't remember more than present atoms ! mem = numbercarbon else mem = memory end if index = 0 error = 0 localradius = 0.0 nonlocalradius = 0.0 smallestradius = 0.0 allocate(coplanarity(4,4)) allocate(N(4,4)) allocate(Z(5,5)) if(thickness .ne. 0.0) then ! Statistical classification of the radii possible now write(*,'(" Second run. Statistics running...")') write(*, *) write(*, *) distribution = 0.0 end if ! Checking all possible combination of 4 points do i = 1 , numbercarbon - 3 do j = i + 1 , numbercarbon - 2 do k = j + 1 , numbercarbon - 1 do l = k + 1, numbercarbon ! coplanarity check -> (d-a).((b-a)\times(c-a)) = 0 coplanarity = reshape((/alphacarbons(i)%coordinates, 1.0, & alphacarbons(j)%coordinates, 1.0, & alphacarbons(k)%coordinates, 1.0, & alphacarbons(l)%coordinates, 1.0/), (/4, 4/)) ! radii calculation D(1:4,1:4) = 0.0 D(1,2) = (sum((alphacarbons(i)%coordinates - & alphacarbons(j)%coordinates)**2)) D(1,3) = (sum((alphacarbons(i)%coordinates - & alphacarbons(k)%coordinates)**2)) D(1,4) = (sum((alphacarbons(i)%coordinates - & alphacarbons(l)%coordinates)**2)) D(2,3) = (sum((alphacarbons(j)%coordinates - & alphacarbons(k)%coordinates)**2)) D(2,4) = (sum((alphacarbons(j)%coordinates - & alphacarbons(l)%coordinates)**2)) D(3,4) = (sum((alphacarbons(k)%coordinates - & alphacarbons(l)%coordinates)**2)) N(1:4,1:4) = 0.0 Z(1:5,1:5) = 0.0 N = D + transpose(D) Z(1:4,1:4) = N Z(1:4,5) = 1.0 Z(5,1:4) = 1.0 detn = determinant(N) detz = determinant(Z) if (detz == 0) then ! The points are coplanar rquadrat = -1 else rquadrat = (- detn / (2 * detz)) end if if ((abs(determinant(coplanarity)) == 0) .OR. (rquadrat <0)) then ! This may also occur if the points are almost coplanar, so that the radius of ! the sphere is quite large if (talkmode) then write(*,'(" Problem with alpha carbons : ", i10, ", ", i10, ", ", i10, ", ", i10)') i, j, k, l write(*,'(" Determinant of Z equals: ", f10.3)') detz write(*,'(" Coplanarity check determinant equals: ", f10.3)') determinant(coplanarity) write(*,*) end if r3 = 0.0 error = error + 1 else r3 = sqrt(rquadrat) if(thickness .ne. 0.0) then ! This is the second run : ! Classification of the found radius pos = floor(((r3/thickness) - 1) * 100) + 1 if (pos .gt. limit) then ! Found radius is bigger than the given limit pos = limit end if distribution(pos) = distribution(pos) + 1 end if combination(1) = i ! This is due to Fortran 90 compatibility combination(2) = j combination(3) = k combination(4) = l do m = 1, 4 ! Completing for all four points s = combination(m) do q = 1, mem if (r3 .lt. smallestradius(s,q)) then ! A smaller radius has been found do r = mem, q+1, -1 ! Shifting the radius vector for insertion of the ! new radius smallestradius(s,r) = smallestradius(s,r-1) smallestlocation(s,r,1) = smallestlocation(s,r-1,1) smallestlocation(s,r,2) = smallestlocation(s,r-1,2) smallestlocation(s,r,3) = smallestlocation(s,r-1,3) smallestlocation(s,r,4) = smallestlocation(s,r-1,4) end do ! Saving the new radius smallestradius(s,q) = r3 smallestlocation(s,q,1) = i smallestlocation(s,q,2) = j smallestlocation(s,q,3) = k smallestlocation(s,q,4) = l exit end if if (smallestradius(s,q) == 0.0) then ! This is the first radius for the given point smallestradius(s,q) = r3 smallestlocation(s,q,1) = i smallestlocation(s,q,2) = j smallestlocation(s,q,3) = k smallestlocation(s,q,4) = l exit end if end do end do end if if (((l-k) + (k-j) + (j-i)) == 3) then ! The found radius is local if ((localradius(j) > r3) .or. (localradius(j) == 0.0)) then localradius(j) = r3 end if ! The local radius is attributed to the both inner points localradius(k) = r3 else ! The radius is not an local one combination(1) = i combination(2) = j combination(3) = k combination(4) = l do m = 1, 4 o = combination(m) if ((nonlocalradius(o) == 0) .or. (nonlocalradius(o) > r3)) then nonlocalradius(o) = r3 nonlocallocation(o,1) = i nonlocallocation(o,2) = j nonlocallocation(o,3) = k nonlocallocation(o,4) = l end if end do end if index = index + 1 end do end do end do end do if (talkmode) then write (*, '(" ", i10, " spheres have been calculated.")') index write (*, '(" For ", i10, " spheres the determination of its radius was not possible.")') error write (*, *) end if deallocate(coplanarity) deallocate(N) deallocate(Z) allocate (globalradius(numbercarbon)) allocate (globallocation(numbercarbon,4)) ! Finalizing do i = 1, numbercarbon ! Which class of radii provides the smallest radius? radii(1) = localradius(i) radii(2) = nonlocalradius(i) selector = minloc(radii, mask = radii .ne. 0.0) location = selector(1) ! Global radius lists the smallest circle for each point globalradius(i) = radii(location) select case (location) ! Indexes of the points realizing the smallest circle case (1) globallocation(i,1) = i - 1 globallocation(i,2) = i globallocation(i,3) = i + 1 globallocation(i,4) = i + 2 case (2) globallocation(i,:) = nonlocallocation(i,:) end select end do allocate (f(numbercarbon)) allocate (g(numbercarbon)) f = 0.0 g = 0.0 do i = 2, numbercarbon - 2 ! Calculation of the f and g factor f(i) = nonlocalradius(i)/localradius(i) g(i) = globalradius(i)/localradius(i) end do deallocate (nonlocalradius, nonlocallocation) deallocate (firstlocal) end subroutine spheres end module distance ! ********************************************************************** module loading implicit none contains subroutine import (inputfile, alphacarbons, atomtype, resolution) ! Loads the preprocessed input data files use var implicit none real :: resolution integer :: i, io, numbercarbons, temp character (LEN=18), intent(in) :: inputfile ! Compatibility with Fortran 90 character (LEN=4) :: atomtype character (LEN=12) :: dummy character (LEN=13) :: dummy2 character (LEN=17) :: dummy3 character (LEN=35) :: test type(carbons), pointer, dimension(:) :: alphacarbons OPEN (UNIT=2, FILE=inputfile, FORM='FORMATTED', STATUS='OLD') ! Reading of the header line and the resolution value rewind 2 read (2, '(a35)') test if (test .ne. 'Atom cartesian coordinates from PDB') then write(*, '("This is not an atom coordinates file")') return end if read (2, *) read (2, *) read (2, '(a13, a3)') dummy2, atomtype read (2, '(a17, f4.0)') dummy3, resolution read (2, *) read (2, *) read (2, *) numbercarbons = 0 do ! Counting the numbers of atoms read (2, '(a12)', IOSTAT = io) dummy if (io == 0) then numbercarbons = numbercarbons + 1 else exit end if end do allocate (alphacarbons(numbercarbons)) rewind 2 do i = 1, 8 read (2, *) end do do i = 1, numbercarbons ! Reading the coordinates read (2, '(i5, 3f10.3, i9, i6, a1)' ) temp, & alphacarbons(i)%coordinates(1), & alphacarbons(i)%coordinates(2), & alphacarbons(i)%coordinates(3), & alphacarbons(i)%residue, temp, & alphacarbons(i)%chainident end do CLOSE (UNIT=2, STATUS='KEEP') end subroutine import end module loading ! ********************************************************************** program calculator ! Main program use loading use var use distance use parameters implicit none character (LEN=12) :: export, dummy character (LEN=7) :: exportfolder, machinefolder character (LEN=5) :: subfolder, humanfolder character (LEN=6) :: importfolder character (LEN=12), allocatable, dimension(:) :: filename character, dimension(:), pointer :: name character (LEN=4) :: extension, atomtype real :: resolution, thickness real(8) :: initial, final real, pointer, dimension(:) :: f, g, localradius real, pointer, dimension(:) :: globalradius real, pointer, dimension(:,:) :: alphacircles, smallestradius type(carbons), pointer, dimension(:) :: alphacarbons type(carbons), allocatable, dimension(:) :: allalphacarbons, oldcarbons integer :: i, j, k, numberfiles, io integer :: position, length, dislocation integer :: a, b, c, d, what, sourceflag integer :: fitflag, numbercarbon, mem integer, allocatable, dimension(:,:) :: comb integer, pointer, dimension(:,:) :: globallocation integer, pointer, dimension(:,:,:) :: smallestlocation integer(4), dimension(limit) :: distribution logical, dimension(:), allocatable :: goodfiles exportfolder = 'results' ! name of the subdirectory where the results are stored to ! Adjust the above field length of the string !!! importfolder = 'source' ! name of the subfolder for ASCII export files ! Adjust the above field length of the string !!! humanfolder = 'human' machinefolder = 'machine' if (talkmode) then write (*, *) write (*, '(" Computational physics project")') write (*, '(" Christian Schlatter, student, physics, 4th year")') write (*, *) releaseinfo write (*, '(" ****************************************************************")') write (*, *) end if ! sourceflag = 1 ! > TEMPORARY ! allocate (filename(1)) ! > TEMPORARY ! filename(1) = '1acc.hel' ! > TEMPORARY ! numberfiles = 1 ! > TEMPORARY write (*, '(" Global curvature analyser")') write (*, *) write (*, '(" How many proteins would you like to process ?")') write (*, *) write (*, '(" -> Type 1 for a single input file or")') write (*, '(" -> Type 2 for a list of files provided in Source.ini")') write (*, *) write (*, '(" Please choose your option now : ")', advance = 'no') sourceflag = 2 ! < TEMPORARY ! read (*,*) sourceflag ! > TEMPORARY write (*, *) if (sourceflag == 2) then ! Batch mode using Source.ini write (*, '(" Opening Source.ini...")') write (*, *) open (UNIT=1, FILE='Source.ini', FORM='FORMATTED', STATUS='OLD') rewind 1 numberfiles = 0 do ! Determination of the number of input files read (1, '(a12)', IOSTAT = io) dummy if (io == 0) then numberfiles = numberfiles + 1 else exit end if end do allocate (filename(numberfiles)) rewind 1 do i = 1, numberfiles ! Determination of the input file names read (1, '(a12)') filename(i) write (*,*) filename(i) end do close (UNIT=1, STATUS='KEEP') write (*,*) write (*, '(" You have chosen ", i3, " files to process")') numberfiles write (*,*) else allocate (filename(1)) ! User provided input file name write (*, '(" Please provide the filename of the import file: ")', advance = 'no') read (*, *) filename(1) write (*, *) numberfiles = 1 endif write (*, '(" What would you like to fit into the protein ?")') write (*, *) write (*, '(" -> Type 1 for fitting circles only")') write (*, '(" -> Type 2 for fitting both - circles and spheres")') write (*, *) write (*, '(" Please choose your option now : ")', advance = 'no') fitflag = 2 ! < TEMPORARY ! read (*,*) fitflag ! > TEMPORARY write (*, *) allocate (goodfiles(numberfiles)) goodfiles = .true. do i = 1, numberfiles if (talkmode) then write (*, *) write (*, '(" Extracting structure information from file ", a12)') filename(i) write (*, '(" *****************************************************************")') write (*, *) end if write (*,'(" Processing now file ", a12)') filename(i) write (*, *) nullify(alphacarbons) if (associated(alphacarbons)) then deallocate(alphacarbons) end if ! Reading the input file call import(importfolder//'/'//filename(i), alphacarbons, atomtype, resolution) numbercarbon = size(alphacarbons) if (memory .gt. numbercarbon) then mem = numbercarbon else mem = memory end if allocate (smallestradius(numbercarbon, mem)) allocate (smallestlocation(numbercarbon,mem,3)) if (talkmode) then write (*, '(" Fitting circles through the alpha carbons")') write (*, '(" *****************************************************************")') write (*, *) end if thickness = 0.0 ! Fitting circles through three atoms call circles (alphacarbons, f, g, localradius, globalradius, & globallocation, smallestradius, smallestlocation, & distribution, thickness) ! Circular thickness of the macromolecule thickness = minval(globalradius) if (secondrun) then ! Doing a second run for statistics on the smallest circles deallocate(localradius, globalradius, f, g, globallocation) call circles (alphacarbons, f, g, localradius, globalradius, & globallocation, smallestradius, smallestlocation, & distribution, thickness) end if position = scan(filename(i),'.') if (filename(i)(position+1:position+1) == 'h') then ! An alpha helix input file has been used subfolder = 'helix' else ! An backbone input file has been used subfolder = 'whole' end if ! Saving the circles to disk filename(i)(position+1:position+1) = 'c' filename(i)(position+2:position+2) = 'i' filename(i)(position+3:position+3) = 'r' open (UNIT=3, FILE=exportfolder//'/'//subfolder//'/'//humanfolder//'/'//filename(i), FORM='FORMATTED', STATUS='REPLACE') open (UNIT=4, FILE=exportfolder//'/'//subfolder//'/'//machinefolder//'/'//filename(i), FORM='FORMATTED', STATUS='REPLACE') ! Making the output file header write(3, '("Circling of proteins")') write(4, '("Circles")') write(3, '("C-alpha-circle project. File created using Calculator.")') write(3, '(a25)', advance='no') releaseinfo write(4, '(a25)') releaseinfo write(3, '(" © 2002 EPFL Lausanne")') write(3, '("Atom type : ", a4)') atomtype write(4, '(a4)') atomtype write(3, '("Resolution [Å] : ", f4.2)') resolution write(4, '(f4.2)') resolution write(3, '("Circular thickness [Å] : ", f4.2)') thickness write(4, '(f4.2)') thickness write(3, '("Thickness combination : ", i3, ", ", i3, " and ", i3)') & globallocation(minloc(globalradius),1), & globallocation(minloc(globalradius),2), & globallocation(minloc(globalradius),3) write(4, '(3i4)') globallocation(minloc(globalradius),1), & globallocation(minloc(globalradius),2), & globallocation(minloc(globalradius),3) write(3, '("Memory : ", i3)') mem write(4, '(i3)') mem if (secondrun) then ! Saving statistics results to disk filename(i)(position+1:position+1) = 'c' filename(i)(position+2:position+2) = 'd' filename(i)(position+3:position+3) = 't' open (UNIT=5, FILE=exportfolder//'/'//subfolder//'/'//machinefolder//'/'//filename(i), FORM='FORMATTED', STATUS='REPLACE') do j = 1, limit ! Distribution of the smallest radii write(5, '(f10.2, i10)') ((j / 100.0) + 1.0), distribution(j) end do close (UNIT=5, STATUS='KEEP') end if write(3, *) write(4, *) write(3, '("Entry Local f-ratio Dis")', advance = 'no') do j = 1, mem write(3, '(" ", i2)', advance = 'no') j write(3, '(". global radius ")', advance = 'no') end do write(3, '(" Residue")') write(3, *) ! Writing data to file do j = 1, (size (alphacarbons)) a = minval(globallocation(j,:)) c = maxval(globallocation(j,:)) b = minval(globallocation(j,:), MASK = globallocation(j,:) .NE. a) ! Spacing of the indexes dislocation = (c - b) + (b - a) write (3, '(i5, 2f10.3, i5)', advance = 'no') j, & localradius(j)/thickness, f(j), dislocation write (4, '(i5, 2f10.3, i5)', advance = 'no') j, & localradius(j)/thickness, f(j), dislocation do k = 1, mem write(3, '(" (", f5.3, ";", 3i4, ") ")', advance = 'no') & smallestradius(j,k)/thickness, & smallestlocation(j,k,1), & smallestlocation(j,k,2), & smallestlocation(j,k,3) write(4, '(f7.3)', advance = 'no') smallestradius(j,k)/thickness end do write (3, '(i8)' ) alphacarbons(j)%residue write (4, '(i8)' ) alphacarbons(j)%residue end do write (*, '(" -> This protein has a thickness of circularily ", f5.3, " Angstroem.")') thickness write(*,*) close (UNIT=3, STATUS='KEEP') close (UNIT=4, STATUS='KEEP') deallocate(localradius, globalradius, f, g, globallocation, smallestradius) deallocate(smallestlocation) if (fitflag == 2) then ! Fitting spheres through 4 atoms of the macromolecule allocate (smallestradius(numbercarbon, mem)) allocate (smallestlocation(numbercarbon,mem,4)) if (talkmode) then write (*, '(" Fitting spheres through the alpha carbons")') write (*, '(" *****************************************************************")') write (*, *) end if thickness = 0.0 call spheres (alphacarbons, f, g, localradius, globalradius, & globallocation, smallestradius, smallestlocation, & distribution, thickness) ! Spherical thickness of the macromolecule write(*, '("First run successfully completed...")') thickness = minval(globalradius) if (secondrun) then ! Doing a second run for statistics on the smallest spheres deallocate(localradius, globalradius, f, g, globallocation) call spheres (alphacarbons, f, g, localradius, globalradius, & globallocation, smallestradius, smallestlocation, & distribution, thickness) end if position = scan(filename(i),'.') ! Saving the spheres to disk filename(i)(position+1:position+1) = 's' filename(i)(position+2:position+2) = 'p' filename(i)(position+3:position+3) = 'h' open (UNIT=3, FILE=exportfolder//'/'//subfolder//'/'//humanfolder//'/'//filename(i), FORM='FORMATTED', STATUS='REPLACE') open (UNIT=4, FILE=exportfolder//'/'//subfolder//'/'//machinefolder//'/'//filename(i), FORM='FORMATTED', STATUS='REPLACE') ! Making the output file header write(3, '("Sphering of proteins")') write(4, '("Spheres")') write(3, '("C-alpha-circle project. File created using Calculator.")') write(3, '(a25)', advance='no') releaseinfo write(4, '(a25)') releaseinfo write(3, '(". © 2002 EPFL Lausanne")') write(3, '("Atom type : ", a4)') atomtype write(4, '(a4)') atomtype write(3, '("Resolution [Å] : ", f4.2)') resolution write(4, '(f4.2)') resolution write(3, '("Spherical thickness [Å] : ", f4.2)') thickness write(4, '(f4.2)') thickness write(3, '("Thickness combination : ", i3, ", ", i3, ", ", i3, " and ", & i3)') globallocation(minloc(globalradius),1), & globallocation(minloc(globalradius),2), & globallocation(minloc(globalradius),3), & globallocation(minloc(globalradius),4) write(4, '(4i4)') globallocation(minloc(globalradius),1), & globallocation(minloc(globalradius),2), & globallocation(minloc(globalradius),3), & globallocation(minloc(globalradius),4) write(3, '("Memory : ", i3)') mem write(4, '(i3)') mem if (secondrun) then ! Saving statistics results to disk filename(i)(position+1:position+1) = 's' filename(i)(position+2:position+2) = 'd' filename(i)(position+3:position+3) = 't' open (UNIT=5, FILE=exportfolder//'/'//subfolder//'/'//machinefolder//'/'//filename(i), FORM='FORMATTED', STATUS='REPLACE') do j = 1, limit ! Distribution of the smallest radii write(5, '(f10.2, i10)') ((j / 100.0) + 1.0), distribution(j) end do close (UNIT=5, STATUS='KEEP') end if write(3, *) write(4, *) write(3, '("Entry Local f-ratio Dis")', advance = 'no') do j = 1, mem write(3, '(" ", i2)', advance = 'no') j write(3, '(". global radius ")', advance = 'no') end do write(3, '(" Residue")') write(3, *) ! Writing data to file do j = 1, (size (alphacarbons)) a = minval(globallocation(j,:)) d = maxval(globallocation(j,:)) b = minval(globallocation(j,:), MASK = globallocation(j,:) .NE. a) c = maxval(globallocation(j,:), MASK = globallocation(j,:) .NE. d) ! Spacing of the indexes dislocation = (d - c) + (c - b) + (b - a) write (3, '(i5, 2f10.3, i5)', advance = 'no') j, & localradius(j)/thickness, f(j), dislocation write (4, '(i5, 2f10.3, i5)', advance = 'no') j, & localradius(j)/thickness, f(j), dislocation do k = 1, mem write(3, '(" (", f5.3, ";", 4i4, ") ")', advance = 'no') & smallestradius(j,k)/thickness, & smallestlocation(j,k,1), & smallestlocation(j,k,2), & smallestlocation(j,k,3), & smallestlocation(j,k,4) write(4, '(f7.3)', advance = 'no') smallestradius(j,k)/thickness end do write (3, '(i8)' ) alphacarbons(j)%residue write (4, '(i8)' ) alphacarbons(j)%residue end do write (*, '(" -> This protein has a thickness of spherically ", f5.3, " Angstroem.")') thickness write(*,*) close (UNIT=3, STATUS='KEEP') close (UNIT=4, STATUS='KEEP') deallocate(localradius, globalradius, f, g, globallocation, smallestradius) deallocate(smallestlocation) end if end do deallocate(filename) deallocate(goodfiles) write(*, '("Program execution finished")') write(*, *) end program calculator