Changeset bd2278d for rmsdfun.f


Ignore:
Timestamp:
09/05/08 11:49:42 (16 years ago)
Author:
baerbaer <baerbaer@…>
Branches:
master
Children:
fafe4d6
Parents:
2ebb8b6
Message:

Reformatting comments and continuation marks.

Fortran 90 and higher use ! to mark comments no matter where they are in the
code. The only valid continuation marker is &.
I also added the SMMP.kdevelop.filelist to the repository to make it easier
to use kdevelop.

git-svn-id: ​svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@12 26dc1dd8-5c4e-0410-9ffe-d298b4865968

File:
1 edited

Legend:

Unmodified
Added
Removed

  • rmsdfun.f

    r2ebb8b6 rbd2278d  
    1 c **************************************************************
    2 c
    3 c This file contains the subroutines: rmsdfun,rmsdopt,fitmol,
    4 c                                     jacobi,rmsinit
    5 c
    6 c Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    7 c                      Shura Hayryan, Chin-Ku
    8 c Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    9 c                      Jan H. Meinke, Sandipan Mohanty
    10 c
    11 c **************************************************************
     1! **************************************************************
     2!
     3! This file contains the subroutines: rmsdfun,rmsdopt,fitmol,
     4!                                     jacobi,rmsinit
     5!
     6! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
     7!                      Shura Hayryan, Chin-Ku
     8! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
     9!                      Jan H. Meinke, Sandipan Mohanty
     10!
     11! **************************************************************
    1212
    1313      real*8 function rmsdfun(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl)
    14 C
    15 C --------------------------------------------------------------
    16 c Wrapping function for calculating rmsd
    17 c
    18 c LIMITATION: requires call of rmsinit BEFORE calling this function
    19 c
    20 c CALLS: rmsdopt
    21 c
    22 c ---------------------------------------------------------------
    23 c
     14!
     15! --------------------------------------------------------------
     16! Wrapping function for calculating rmsd
     17!
     18! LIMITATION: requires call of rmsinit BEFORE calling this function
     19!
     20! CALLS: rmsdopt
     21!
     22! ---------------------------------------------------------------
     23!
    2424      include 'INCL.H'
    2525      include 'INCP.H'
    26 c
    27 c Input
     26!
     27! Input
    2828      dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
    29 c Local
     29! Local
    3030      dimension rm(3,3),av1(3),av2(3)
    3131      call rmsdopt(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl,rm,av1,av2,rssd)
    … …  
    3737      end
    3838
    39 c*******************************************************************
     39!*******************************************************************
    4040      subroutine rmsdopt(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl,
    41      #                   rm,av1,av2,rmsd)
    42 
    43 c ---------------------------------------------------------------
    44 c PURPOSE: root mean square deviation (rmsd) between current SMMP
    45 c          structure and reference atom coordinates 'x,y,zrf()'
    46 c          for range of SMMP residues [ir1,ir2] in molecule 'nml'
    47 c
    48 c          ixat(i) - points to the atom in ref. coords., which is
    49 c                    equivalent to atom i of SMMP structure
    50 c                    (=0 if no equivalent in ref. structure exists)
    51 c
    52 c          isl = 0  : select all heavy atoms
    53 c          isl = 1  : backbone atoms n,ca,c
    54 c          isl = 2  : only ca atoms
    55 c
    56 c CALLS:   fitmol [S.K.Kearsley, Acta Cryst. 1989, A45, 208-210]
    57 c
    58 c      NB  uncomment last lines in 'fitmol' to return coordinates
    59 c          in 'x2' after fitting the ref. str. onto SMMP structure
    60 c ----------------------------------------------------------------
     41     &                   rm,av1,av2,rmsd)
     42
     43! ---------------------------------------------------------------
     44! PURPOSE: root mean square deviation (rmsd) between current SMMP
     45!          structure and reference atom coordinates 'x,y,zrf()'
     46!          for range of SMMP residues [ir1,ir2] in molecule 'nml'
     47!
     48!          ixat(i) - points to the atom in ref. coords., which is
     49!                    equivalent to atom i of SMMP structure
     50!                    (=0 if no equivalent in ref. structure exists)
     51!
     52!          isl = 0  : select all heavy atoms
     53!          isl = 1  : backbone atoms n,ca,c
     54!          isl = 2  : only ca atoms
     55!
     56! CALLS:   fitmol [S.K.Kearsley, Acta Cryst. 1989, A45, 208-210]
     57!
     58!      NB  uncomment last lines in 'fitmol' to return coordinates
     59!          in 'x2' after fitting the ref. str. onto SMMP structure
     60! ----------------------------------------------------------------
    6161
    6262      include 'INCL.H'
    6363      include 'INCP.H'
    6464
    65 c-------------------------------------------------------- input
     65!-------------------------------------------------------- input
    6666      dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
    67 c-------------------------------------------------------- output
     67!-------------------------------------------------------- output
    6868      dimension rm(3,3),av1(3),av2(3)
    69 c-------------------------------------------------------- local
     69!-------------------------------------------------------- local
    7070      dimension x1(3,mxat),x2(3,mxat)
    7171      character*4 atnm
    … …  
    9999                if ( isl.eq.0
    100100
    101      #               .or.
    102 
    103      #              (isl.eq.1.and.(index(atnm,'n   ').gt.0 .or.
    104      #                             index(atnm,'ca  ').gt.0 .or.
    105      #                             index(atnm,'c   ').gt.0 ))
    106      #               .or.
    107 
    108      #              (isl.eq.2.and.index(atnm,'ca  ').gt.0)
    109 
    110      #          ) then
     101     &               .or.
     102
     103     &              (isl.eq.1.and.(index(atnm,'n   ').gt.0 .or.
     104     &                             index(atnm,'ca  ').gt.0 .or.
     105     &                             index(atnm,'c   ').gt.0 ))
     106     &               .or.
     107
     108     &              (isl.eq.2.and.index(atnm,'ca  ').gt.0)
     109
     110     &          ) then
    111111
    112112                  n=n+1
    … …  
    136136      return
    137137      end
    138 c *********************************************
     138! *********************************************
    139139      subroutine fitmol(n,x1,x2, rm,a1,a2,rmsd)
    140 c      real*8 function fitmol(n,x1,x2)
    141 
    142 c .......................................................
    143 c PURPOSE: compute RMSD of n positions in x1(3,) & x2(3,)
    144 c          [S.K.Kearsley Acta Cryst. 1989,A45,208-210]
    145 c
    146 c CALLS: jacobi
    147 c .......................................................
    148 cf2py intent(out) rmsd
     140!      real*8 function fitmol(n,x1,x2)
     141
     142! .......................................................
     143! PURPOSE: compute RMSD of n positions in x1(3,) & x2(3,)
     144!          [S.K.Kearsley Acta Cryst. 1989,A45,208-210]
     145!
     146! CALLS: jacobi
     147! .......................................................
     148!f2py intent(out) rmsd
    149149                   
    150150      include 'INCL.H'
    151 c      implicit real*8 (a-h,o-z)
    152 c      implicit integer*4 (i-n)
     151!      implicit real*8 (a-h,o-z)
     152!      implicit integer*4 (i-n)
    153153 
    154 c ------------------------------------------- input/output
     154! ------------------------------------------- input/output
    155155      dimension x1(3,mxat),x2(3,mxat)
    156 c -------------------------------------------------- local
     156! -------------------------------------------------- local
    157157      dimension e(4),q(4,4),v(4,4),dm(3),dp(3),a1(3),a2(3),rm(3,3)
    158158
    159159      dn=dble(n)
    160 c ------------------- average of coordinates
     160! ------------------- average of coordinates
    161161      do i=1,3
    162162        a1(i) = 0.d0
    … …  
    169169        a2(i) = a2(i)/dn
    170170      enddo
    171 c ------------------------- compile quaternion
     171! ------------------------- compile quaternion
    172172      do i=1,4
    173173        do j=1,4
    … …  
    208208        enddo
    209209      enddo
    210 c ------------------------------ eigenvalues & -vectors
     210! ------------------------------ eigenvalues & -vectors
    211211      ndim4=4
    212212      call jacobi(q,ndim4,e,v)
    213 c --------------------------- lowest eigenvalue
     213! --------------------------- lowest eigenvalue
    214214      im=1
    215215      em=e(1)
    … …  
    223223      rmsd = sqrt(em/dn)
    224224
    225 c ================= uncomment following lines to fit molecule 2 onto 1
    226 
    227 c ---------------------------------------------------rotation matrix
     225! ================= uncomment following lines to fit molecule 2 onto 1
     226
     227! ---------------------------------------------------rotation matrix
    228228      rm(1,1) = v(1,im)**2+v(2,im)**2-v(3,im)**2-v(4,im)**2
    229229      rm(1,2) = 2.d0*( v(2,im)*v(3,im)-v(1,im)*v(4,im) )
    … …  
    236236      rm(3,3) = v(1,im)**2+v(4,im)**2-v(2,im)**2-v(3,im)**2
    237237
    238 c      do i=1,n
    239 c        do j=1,3
    240 c          dm(j) = x2(j,i) - a2(j)
    241 c        enddo
    242 c        do j=1,3
    243 c          dp(j) = a1(j)
    244 c          do k=1,3
    245 c            dp(j) = dp(j) + rm(j,k) * dm(k)
    246 c          enddo
    247 c          x2(j,i) = dp(j)
    248 c        enddo
    249 c      enddo
    250 
    251 c      fitmol=rmsd
     238!      do i=1,n
     239!        do j=1,3
     240!          dm(j) = x2(j,i) - a2(j)
     241!        enddo
     242!        do j=1,3
     243!          dp(j) = a1(j)
     244!          do k=1,3
     245!            dp(j) = dp(j) + rm(j,k) * dm(k)
     246!          enddo
     247!          x2(j,i) = dp(j)
     248!        enddo
     249!      enddo
     250
     251!      fitmol=rmsd
    252252
    253253      return
    254254      end
    255 c ******************************
     255! ******************************
    256256      subroutine jacobi(a,n,d,v)
    257257
    258 c ......................................................
    259 c PURPOSE: for given symmetric matrix 'a(n,n)
    260 c          compute eigenvalues 'd' & eigenvectors 'v(,)'
    261 c
    262 c  [W.H.Press,S.A.Teukolsky,W.T.Vetterling,
    263 c   B.P.Flannery, Numerical Recipes in FORTRAN,
    264 c   Cambridge Univ. Press, 2nd Ed. 1992, 456-462]
    265 c
    266 c CALLS: none
    267 c
    268 c ......................................................
    269 
    270 cf2py intent(out) d
    271 cf2py intent(out) v
     258! ......................................................
     259! PURPOSE: for given symmetric matrix 'a(n,n)
     260!          compute eigenvalues 'd' & eigenvectors 'v(,)'
     261!
     262!  [W.H.Press,S.A.Teukolsky,W.T.Vetterling,
     263!   B.P.Flannery, Numerical Recipes in FORTRAN,
     264!   Cambridge Univ. Press, 2nd Ed. 1992, 456-462]
     265!
     266! CALLS: none
     267!
     268! ......................................................
     269
     270!f2py intent(out) d
     271!f2py intent(out) v
    272272      parameter (NMAX=500)
    273273     
    … …  
    276276
    277277      real*8 a(n,n),d(n),v(n,n),
    278      #       c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX),smeps
     278     &       c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX),smeps
    279279
    280280      smeps=1.0d-6
    … …  
    318318            if((i.gt.4).and.(abs(d(ip))+
    319319
    320      #g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then
     320     &g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then
    321321              a(ip,iq)=0.d0
    322322
    … …  
    393393      end
    394394
    395 c ***********************************************************
     395! ***********************************************************
    396396
    397397      subroutine rmsinit(nml,string)
    398 c
    399 c------------------------------------------------------------------------------
    400 c Reads in pdb-file 'string' into INCP.H and initalizes
    401 c the files that 'rmdsopt' needs to calculate the rmsd
    402 c of a configuration with the pdb-configuration
    403 C
    404 c CALLS: pdbread,atixpdb
    405 c
    406 c ----------------------------------------------------------------------------
    407 c
     398!
     399!------------------------------------------------------------------------------
     400! Reads in pdb-file 'string' into INCP.H and initalizes
     401! the files that 'rmdsopt' needs to calculate the rmsd
     402! of a configuration with the pdb-configuration
     403!
     404! CALLS: pdbread,atixpdb
     405!
     406! ----------------------------------------------------------------------------
     407!
    408408      include 'INCL.H'
    409409      include 'INCP.H'
    … …  
    412412 
    413413      if(string.eq.'smmp') then
    414 c
    415 c Compare with a smmp-structure
    416 c
     414!
     415! Compare with a smmp-structure
     416!
    417417        do i=iatrs1(irsml1(nml)),iatrs2(irsml2(nml))
    418418         if(nmat(i)(1:1).ne.'h') then
    … …  
    422422         end if
    423423        enddo
    424 c
     424!
    425425       else
    426 c
    427 c Reference structure is read in from pdb-file
    428 c
     426!
     427! Reference structure is read in from pdb-file
     428!
    429429         call pdbread(string,ier)
    430430         if(ier.ne.0) stop
    431431         call atixpdb(nml)
    432 c
     432!
    433433      end if
    434434      print *,'RMSD initialized with ',string
Note: See TracChangeset for help on using the changeset viewer.