Changeset 32289cd for esolan.f

Timestamp:

11/19/09 11:29:41 (14 years ago)

Author:

baerbaer <baerbaer@…>

Branches:

master

Children:

38d77eb

Parents:

6650a56

Message:

Explicitly declare variables.

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

File:

• esolan.f (modified) (39 diffs)

esolan.f

@@ -4 +4 @@
 !
 ! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
-!                      Shura Hayryan, Chin-Ku 
+!                      Shura Hayryan, Chin-Ku
 ! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
 !                      Jan H. Meinke, Sandipan Mohanty
@@ -14 +14 @@
 ! -----------------------------------------------------------------
 !
-! Calculates the solvation energy of the protein with 
-! solavtion parameters model E=\sum\sigma_iA_i. 
+! Calculates the solvation energy of the protein with
+! solavtion parameters model E=\sum\sigma_iA_i.
 ! The solvent accessible surface area per atom
-! and its gradients with respect to the Cartesian coordinates of 
+! and its gradients with respect to the Cartesian coordinates of
 ! the atoms are calculated using the projection method described in
 !
@@ -25 +25 @@
 !
 ! INPUT: nmol - the order number of the protein chain.
-!      The atomic coordinates, radii and solvation parameters  
+!      The atomic coordinates, radii and solvation parameters
 !      are taken from the global arrays of SMMP
 !
-! OUTPUT: global array gradan(mxat,3) contains the gradients of 
+! OUTPUT: global array gradan(mxat,3) contains the gradients of
 !         solvation energy per atom
 !         local array as(mxat) contains accessible surface area per atom
@@ -42 +42 @@
 
       include 'INCL.H'
-      
-      
+
+
       integer nmol, ks0, ks2
 Cf2py intent(in) nmol
-      
+
       parameter (ks0=mxat*(mxat+1))
       parameter (ks2=mxat+mxat)
@@ -53 +53 @@
 
       integer neib, neibp, ivrx, iv, ial, i, ij, j, iii, idi,ivr, ifu,
-     &          ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll, 
+     &          ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll,
      &          nb, numat, nyx
       real*8 vertex, ax, pol, as, ayx, ayx1, probe, dd, ddat, dadx, gp,
      &          dalp,dbet, daalp, dabet, vrx, dv, dx, dy, dz, dt, di,
-     &          dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold, 
+     &          dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold,
      &          energy, sss, a1, a, ab1, a2, aa, ac2, ac1, ab2, am2,
      &          am, aia, ak, alp, am1, ap1, amibe, amabe, amaal, amial,
@@ -63 +63 @@
      &          ba, bcd, bc, bet, c1, c, c2, cc, sfcg, caba, cg, cfsg,
      &          cfcg, daba, ct, cs, d1, d2, dbe, dal, dcotbma, dcbet,
-     &          dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp, 
-     &          dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22, 
+     &          dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp,
+     &          dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22,
      &          r22p, r2, pom, prob, r, r42p, ratom, rr, s, sfsg, sf,
-     &          sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv, 
+     &          sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv,
      &          zz, zv
       dimension neib(0:mxat),neibp(0:mxat),ivrx(ks0)
@@ -78 +78 @@
       integer ta2(0:mxat),ta3(0:mxat),fullarc(0:mxat),al(0:ks2)
       real*8 neibor(mxat,4)
-      
+
       real*8, dimension(:, :, :), allocatable :: grad
-      
+
       allocate(grad(mxat, mxat, 3))
-      
-      
+
+
 !     open(3,file='solvation.dat')! Output file with solvation data
 
@@ -90 +90 @@
       if (nmol.eq.0) then
          numat=iatrs2(irsml2(ntlml))-iatrs1(irsml1(1))+1 ! Number of atoms
-      else 
+      else
          numat=iatrs2(irsml2(nmol))-iatrs1(irsml1(nmol))+1 ! Number of atoms
       endif
@@ -96 +96 @@
       do i=1,numat
         do j=1,3
-         gradan(i,j)=0.0d0 
+         gradan(i,j)=0.0d0
       end do
       end do
@@ -102 +102 @@
       do i=1,numat
          xold(i)=xat(i)! Redefine the coordinates just for safety
-         yold(i)=yat(i)! 
-         zold(i)=zat(i)! 
-      pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main' 
+         yold(i)=yat(i)!
+         zold(i)=zat(i)!
+      pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main'
       As(i)=pi4*pol(i)! Initially the whole surface of the atom is accessible.
       end do
-      
+
              do iii=1,numat
               do jjj=1,numat
@@ -120 +120 @@
 !***************************
 
-      do 1 i=1,numat ! Lop over atoms "i"           
+      do 1 i=1,numat ! Lop over atoms "i"
 ! ----------------------------------------------------
-         
-       R=rvdw(i) 
+
+       R=rvdw(i)
        jj=0
         do j=1,numat !Find the neighbours of "i"
@@ -132 +132 @@
          write(*,*)'ERROR in data: centres of two atoms coincide!'
          write(*,*)i,j,xold(i),yold(i),zold(i),rvdw(i),
-     &                   xold(j),yold(j),zold(j),rvdw(j) 
+     &                   xold(j),yold(j),zold(j),rvdw(j)
          stop 'Centres of atoms coincide!'
        endif
@@ -146 +146 @@
          end if
          end do!Neighbours
-       
+
        neib(0)=jj   ! The number of neighbors of "i"
 
@@ -197 +197 @@
 !
         uga=grnd() ! Random \gamma angle
-        ufi=grnd() ! Random \pi angle 
+        ufi=grnd() ! Random \pi angle
         uga=pi*uga/2
         ufi=pi*ufi/2
@@ -219 +219 @@
         end do
       end do
-!       
-      
+!
+
        if(jjj.ne.0) then ! Rotation is necessary
          sfsg=sf*sg
@@ -243 +243 @@
       do jj=1,neib(0)
       j=neib(jj)
-       neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+      
+       neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+
      &            (ddat(j,4)-R)**2-pol(j)           ! a
         neibor(j,2)=-pom*ddat(j,2)                  ! b
@@ -259 +259 @@
        do while(k.gt.1)               ! B
         k=k-1
-! Analyse mutual disposition of every pair of neighbours                                                      
+! Analyse mutual disposition of every pair of neighbours
         do 13 L=neib(0),k+1,-1        ! A
 
@@ -283 +283 @@
           elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.      ! a01 02
      &    -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
-! The circle neib(k) encloses circle neib(L) and the later is discarded 
+! The circle neib(k) encloses circle neib(L) and the later is discarded
            neib(0)=neib(0)-1
            do j=L,neib(0)
@@ -326 +326 @@
           elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.  ! a03 02
      &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
-! Don't exclude neib(k) 
+! Don't exclude neib(k)
            neib(0)=neib(0)-1
            do j=k,neib(0)
@@ -370 +370 @@
           elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.   ! a07 02
      &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
-! discard the circle neib(L) 
+! discard the circle neib(L)
            neib(0)=neib(0)-1
            do j=L,neib(0)
@@ -385 +385 @@
            p2=(b1-b2)*pom
 !-- Assign t and s coordinates and the order number of circles
-           vertex(iv,1)=(t+p1)/am 
+           vertex(iv,1)=(t+p1)/am
            vertex(iv,2)=(s-p2)/am
            vertex(iv,3)=neib(k)
@@ -394 +394 @@
            vertex(iv,3)=neib(k)
            vertex(iv,4)=neib(L)
-!--  
+!--
           endif                                                 ! 10
 
@@ -491 +491 @@
             dbe=amabe-amibe
             if(dal.lt.pi) then
-             As(i)=R22*(pi-dal)  
+             As(i)=R22*(pi-dal)
             elseif(dbe.lt.pi) then
              As(i)=R22*(pi-dbe)
@@ -510 +510 @@
        enddo
        do j=1,neib(0)
-! "iv" is the number of intersection point. 
+! "iv" is the number of intersection point.
         do k=1,iv
        if(neibp(j).eq.vertex(k,3).or.neibp(j).eq.vertex(k,4)) goto 30
@@ -532 +532 @@
        do k=1,ifu
          a=neibor(fullarc(k),1)
-         if(a.lt.0.d0) then 
+         if(a.lt.0.d0) then
           aia=-1.d0
          else
@@ -579 +579 @@
          ak=a*a
          daba=dabs(a)
-         if(a.lt.0d0) then 
+         if(a.lt.0d0) then
           aia=-1d0
           else
@@ -679 +679 @@
           do j=1,nyx
 
-!  
+!
 !  Escape 'bad' intersections.
            if(dabs(ayx1(j+1)-ayx1(j)).lt.1d-8) goto 40
@@ -686 +686 @@
            ap1=ct+rr*dcos(prob) ! The middle point of the arc.
            ap2=cs+rr*dsin(prob)
-! Verify if the middle point belongs to a covered part. 
+! Verify if the middle point belongs to a covered part.
 ! If yes then omit.
            do ij=1,ial
@@ -754 +754 @@
              dabet(3)=(a*ac2-ak*c2+wv*daba*a2*dsbet)/am2
              dabet(4)=(2d0*ak*a2)/am2
-               dv(2)=R42*b*vv1      
-               dv(3)=R42*c*vv1      
-               dv(4)=(d-R42*a)*vv1      
-               dv(1)=-dv(4)*R42      
-               dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2      
-             dx(2)=-(d+R42*a)*dv(2)*vv2      
-             dx(3)=-(d+R42*a)*dv(3)*vv2      
-             dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2      
+               dv(2)=R42*b*vv1
+               dv(3)=R42*c*vv1
+               dv(4)=(d-R42*a)*vv1
+               dv(1)=-dv(4)*R42
+               dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2
+             dx(2)=-(d+R42*a)*dv(2)*vv2
+             dx(3)=-(d+R42*a)*dv(3)*vv2
+             dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2
              dy(1)=(-2d0*d+4d0*R42*a)/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
      &             R42*ak)*(aia*vv+daba*dv(1))/ak*vv2
@@ -769 +769 @@
      &             R42*ak)*aia*dv(3)/a*vv2
                dy(4)=-2.d0*a/daba*vv1-(b*b+c*c-2.d0*a*d+2.d0*
-     &             R42*ak)*aia*dv(4)/a*vv2      
+     &             R42*ak)*aia*dv(4)/a*vv2
              dz(1)=-2d0*d/wv/a*vv1-wv*(vv+a*dv(1))/ak*vv2
              dz(2)=b/wv/a*vv1-wv*dv(2)/a*vv2
@@ -820 +820 @@
              dii(4,2)=dii(1,2)*R42
              dii(4,3)=2d0*(ddat(jk,4)+R)*R42
-             
+
              do idi=1,3
               grad(i,jk,idi)=grad(i,jk,idi)+
@@ -826 +826 @@
      &                   di(3)*dii(3,idi)+di(4)*dii(4,idi)
 
-             enddo             
+             enddo
              do idi=1,4
                dta(idi)=5d-1*daalp(idi)*(((yv-zv*((-b*dsbetpalp
@@ -838 +838 @@
              di1(idi)=R2*(aia*daalp(idi)-4d0*xv*dta(idi)/(4d0+tt*tt))
             di2(idi)=R2*(-aia*dabet(idi)-4d0*xv*dtb(idi)/(4d0+tt*tt))
-             enddo             
+             enddo
 
              dii(1,1)=2d0*ddat(kk,2)
@@ -866 +866 @@
 40           continue
           enddo
-         
+
          else ! analyse the case with lines (not necessary after rotation).
 
@@ -920 +920 @@
            enddo
            probe(nyx+1)=ayx1(nyx)+1d0
-           
+
 
            bc=b*b+c*c
@@ -960 +960 @@
 22           continue
             endif
-         enddo       
+         enddo
 11      continue
 
@@ -1017 +1017 @@
 111   continue
 
-!     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',     
+!     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',
 !    & '   gradz   Rad    Atom'
 !     write(3,*)
 
-      
+
       do i=1,numat
         enr=As(i)*sigma(i)
@@ -1038 +1038 @@
 
       esolan=sss
-      
+
       deallocate(grad)
 !     close(3)
@@ -1046 +1046 @@
 200   format(2i4,3f16.6)
 700   format(i5,7f8.3,5x,a4)
-         
+
       end