Changeset b477fe8


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

Fixed handling of nmol = 0 in esolan.

The function energy calls esolan(0) to get the solvent energy if itysol > 0.
Esolan should calculate the solvent energy using all atoms in the system in
this case. I added an if statement at the beginning of esolan() that takes
care of this now. Further rewrote energy() so that it assigns the result
returned by esolan() to teysl, which in turn will be assigned to eysl in the
end.

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

Files:
2 edited

Legend:

Unmodified
Added
Removed

  • energy.f

    r078aff3 rb477fe8  
    7272            esm=esm+enysol(0)
    7373            teysl = teysl+eysl
    74       elseif (itysol.lt.0) then
    75       esm=esm+esolan(0)
     74         else if (itysol.lt.0) then
     75            eysl = esolan(0)
     76            teysl = teysl + eysl
     77            esm = esm + eysl
    7678         else
    7779            eysl=0.d0
    … …  
    143145            esm=esm+enysol(nml)
    144146      elseif (itysol.lt.0) then
    145       esm=esm+esolan(nml)
     147            esm=esm+esolan(nml)
    146148         else
    147149            eysl=0.d0

  • esolan.f

    r078aff3 rb477fe8  
    1 c **************************************************************
    2 c
    3 c This file contains the subroutines: esolan
    4 c
    5 c Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 c                      Shura Hayryan, Chin-Ku
    7 c Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    8 c                      Jan H. Meinke, Sandipan Mohanty
    9 c
    10 c **************************************************************
     1! **************************************************************
     2!
     3! This file contains the subroutines: esolan
     4!
     5! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
     6!                      Shura Hayryan, Chin-Ku
     7! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
     8!                      Jan H. Meinke, Sandipan Mohanty
     9!
     10! **************************************************************
    1111
    1212      real*8 function esolan(nmol)
    1313
    14 c -----------------------------------------------------------------
    15 c
    16 c Calculates the solvation energy of the protein with
    17 c solavtion parameters model E=\sum\sigma_iA_i.
    18 c The solvent accessible surface area per atom
    19 c and its gradients with respect to the Cartesian coordinates of
    20 c the atoms are calculated using the projection method described in
    21 c
    22 c            J Comp Phys 26, 334-343, 2005
    23 c
    24 C USAGE: eysl=esolan(nmol)-Returns the value of the solvation energy
    25 c
    26 c INPUT: nmol - the order number of the protein chain.
    27 c      The atomic coordinates, radii and solvation parameters 
    28 c      are taken from the global arrays of SMMP
    29 c
    30 c OUTPUT: global array gradan(mxat,3) contains the gradients of
    31 c         solvation energy per atom
    32 c         local array as(mxat) contains accessible surface area per atom
    33 c
    34 c Output file "solvation.dat' conatains detailed data.
    35 c
    36 c Correctness of this program was last time rechecked with GETAREA
    37 c
    38 c CALLS: none
    39 c
    40 c ----------------------------------------------------------------------
    41 c TODO Test the solvent energy for multiple molecules
     14! -----------------------------------------------------------------
     15!
     16! Calculates the solvation energy of the protein with
     17! solavtion parameters model E=\sum\sigma_iA_i.
     18! The solvent accessible surface area per atom
     19! and its gradients with respect to the Cartesian coordinates of
     20! the atoms are calculated using the projection method described in
     21!
     22!            J Comp Phys 26, 334-343, 2005
     23!
     24! USAGE: eysl=esolan(nmol)-Returns the value of the solvation energy
     25!
     26! INPUT: nmol - the order number of the protein chain.
     27!      The atomic coordinates, radii and solvation parameters 
     28!      are taken from the global arrays of SMMP
     29!
     30! OUTPUT: global array gradan(mxat,3) contains the gradients of
     31!         solvation energy per atom
     32!         local array as(mxat) contains accessible surface area per atom
     33!
     34! Output file "solvation.dat' conatains detailed data.
     35!
     36! Correctness of this program was last time rechecked with GETAREA
     37!
     38! CALLS: none
     39!
     40! ----------------------------------------------------------------------
     41! TODO Test the solvent energy for multiple molecules
    4242
    4343      include 'INCL.H'
    … …  
    5959     
    6060     
    61 c     open(3,file='solvation.dat')! Output file with solvation data
     61!     open(3,file='solvation.dat')! Output file with solvation data
    6262
    6363      energy=0.0d0 ! Total solvation energy
    6464      sss=0.0d0 ! Total solvent accessible surface area
    65       numat=iatrs2(irsml2(nmol))-iatrs1(irsml1(nmol))+1 ! Number of atoms
     65      if (nmol.eq.0) then
     66         numat=iatrs2(irsml2(ntlml))-iatrs1(irsml1(1))+1 ! Number of atoms
     67      else
     68         numat=iatrs2(irsml2(nmol))-iatrs1(irsml1(nmol))+1 ! Number of atoms
     69      endif
    6670
    6771      do i=1,numat
    … …  
    8791             enddo
    8892
    89 c***************************
    90 c   Start the computations *
    91 c***************************
     93!***************************
     94!   Start the computations *
     95!***************************
    9296
    9397      do 1 i=1,numat ! Lop over atoms "i"           
    94 c ----------------------------------------------------
     98! ----------------------------------------------------
    9599         
    96100       R=rvdw(i)
    … …  
    99103       if(i.ne.j) then
    100104       ddp=(xold(j)-xold(i))**2+(yold(j)-yold(i))**2
    101      #   +(zold(j)-zold(i))**2! Distance between atom i and j
     105     &   +(zold(j)-zold(i))**2! Distance between atom i and j
    102106       if(ddp.lt.1e-10) then
    103107         write(*,*)'ERROR in data: centres of two atoms coincide!'
    104108         write(*,*)i,j,xold(i),yold(i),zold(i),rvdw(i),
    105      #                   xold(j),yold(j),zold(j),rvdw(j)
     109     &                   xold(j),yold(j),zold(j),rvdw(j)
    106110         stop 'Centres of atoms coincide!'
    107111       endif
    … …  
    122126       if(neib(0).eq.0) goto 1 !Finish the atom i if it doesn't have neighbors
    123127
    124 c----
     128!----
    125129       R2=R*R                   ! R square
    126130       R22=R2+R2             ! 2 * R square
    … …  
    140144           ddat(i,4)=0d0
    141145           ddat(i,1)=R
    142 c
    143 c    Verification of the north point of ith sphere
    144 c
     146!
     147!    Verification of the north point of ith sphere
     148!
    145149       jjj=0
    146 c
    147 c      If jjj=0 then - no transformation
    148 c             else the transformation is necessary
    149 c
     150!
     151!      If jjj=0 then - no transformation
     152!             else the transformation is necessary
     153!
    150154       k=neib(1)  ! The order number of the first neighbour
    151 c ddp - the square minimal distance of NP from the neighboring surfaces
     155! ddp - the square minimal distance of NP from the neighboring surfaces
    152156       ddp=dabs((ddat(k,2))**2+(ddat(k,3))**2+
    153      #    (R-ddat(k,4))**2-pol(k))
     157     &    (R-ddat(k,4))**2-pol(k))
    154158
    155159       do j=2,neib(0)
    156160        k=neib(j)
    157161        ddp2=dabs((ddat(k,2))**2+(ddat(k,3))**2+
    158      #      (R-ddat(k,4))**2-pol(k))
     162     &      (R-ddat(k,4))**2-pol(k))
    159163        if(ddp2.lt.ddp) ddp=ddp2
    160164       enddo
    161 c
    162 c   Check whether the NP of ith sphere is too close to the intersection line
    163 c
     165!
     166!   Check whether the NP of ith sphere is too close to the intersection line
     167!
    164168       do while (ddp.lt.1.d-6)
    165169        jjj=jjj+1
    166 c
    167 c       generate grndom numbers
    168 c
     170!
     171!       generate grndom numbers
     172!
    169173        uga=grnd() ! Random \gamma angle
    170174        ufi=grnd() ! Random \pi angle
    … …  
    182186        k=neib(1)
    183187        ddp=dabs((xx-ddat(k,2))**2+(yy-ddat(k,3))**2+
    184      #     (zz-ddat(k,4))**2-pol(k))
     188     &     (zz-ddat(k,4))**2-pol(k))
    185189        do j=2,neib(0)
    186190         k=neib(j)
    187191         ddp2=dabs((xx-ddat(k,2))**2+(yy-ddat(k,3))**2+
    188      #       (zz-ddat(k,4))**2-pol(k))
     192     &       (zz-ddat(k,4))**2-pol(k))
    189193         if(ddp2.lt.ddp) ddp=ddp2
    190194        end do
    191195      end do
    192 c       
     196!       
    193197     
    194198       if(jjj.ne.0) then ! Rotation is necessary
    … …  
    206210       endif
    207211
    208 c  iiiiiiiiiii
     212!  iiiiiiiiiii
    209213
    210214       pom=8.d0*pol(i)
    211215
    212 C In this loop the parameters a,b,c,d for the equation
    213 c      a*(t^2+s^2)+b*t+c*s+d=0 are calculated (see the reference article)
     216! In this loop the parameters a,b,c,d for the equation
     217!      a*(t^2+s^2)+b*t+c*s+d=0 are calculated (see the reference article)
    214218      do jj=1,neib(0)
    215219      j=neib(jj)
    216220       neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+     
    217      #            (ddat(j,4)-R)**2-pol(j)           ! a
     221     &            (ddat(j,4)-R)**2-pol(j)           ! a
    218222        neibor(j,2)=-pom*ddat(j,2)                  ! b
    219223        neibor(j,3)=-pom*ddat(j,3)                  ! c
    220224        neibor(j,4)=4d0*pol(i)*((ddat(j,2))**2+
    221      #      (ddat(j,3))**2+(R+ddat(j,4))**2-pol(j)) ! d
     225     &      (ddat(j,3))**2+(R+ddat(j,4))**2-pol(j)) ! d
    222226       enddo
    223227
    224 c       end of the 1st cleaning
     228!       end of the 1st cleaning
    225229
    226230       iv=0
    … …  
    230234       do while(k.gt.1)               ! B
    231235        k=k-1
    232 C Analyse mutual disposition of every pair of neighbours                                                     
     236! Analyse mutual disposition of every pair of neighbours                                                     
    233237        do 13 L=neib(0),k+1,-1        ! A
    234238
    … …  
    242246          d2=neibor(neib(L),4)/neibor(neib(L),1)
    243247          D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
    244      #      (d1-d2)**2)+(b1*c2-b2*c1)**2
    245 C if D<0 then the circles neib(k) and neib(L) don't intersect
     248     &      (d1-d2)**2)+(b1*c2-b2*c1)**2
     249! if D<0 then the circles neib(k) and neib(L) don't intersect
    246250          if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.      ! a01 01
    247      #     dsqrt(b2*b2+c2*c2-4d0*d2)-dsqrt(b1*b1+c1*c1-4d0*d1)) then ! 04
    248 C Circle neib(L) encloses circle neib(k) and the later is discarded
     251     &     dsqrt(b2*b2+c2*c2-4d0*d2)-dsqrt(b1*b1+c1*c1-4d0*d1)) then ! 04
     252! Circle neib(L) encloses circle neib(k) and the later is discarded
    249253           neib(0)=neib(0)-1
    250254           do j=k,neib(0)
    … …  
    253257           goto 12
    254258          elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.      ! a01 02
    255      #    -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
    256 C The circle neib(k) encloses circle neib(L) and the later is discarded
     259     &    -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
     260! The circle neib(k) encloses circle neib(L) and the later is discarded
    257261           neib(0)=neib(0)-1
    258262           do j=L,neib(0)
    … …  
    260264           enddo
    261265          elseif(D.gt.0d0) then                        ! a01 03
    262 C The circles nieb(L) and neib(k) have two intersection points (IP)
     266! The circles nieb(L) and neib(k) have two intersection points (IP)
    263267           am=2d0*((b1-b2)**2+(c1-c2)**2)
    264268           t=-2d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
    … …  
    280284
    281285         elseif(neibor(neib(k),1).gt.0d0.and.      ! a03
    282      #              neibor(neib(L),1).lt.0d0) then
     286     &              neibor(neib(L),1).lt.0d0) then
    283287          b1=neibor(neib(k),2)/neibor(neib(k),1)
    284288          c1=neibor(neib(k),3)/neibor(neib(k),1)
    … …  
    288292          d2=neibor(neib(L),4)/neibor(neib(L),1)
    289293          D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
    290      #        (d1-d2)**2)+(b1*c2-b2*c1)**2
    291 C if D<0 then neib(k) and neib(L) don't intersect
     294     &        (d1-d2)**2)+(b1*c2-b2*c1)**2
     295! if D<0 then neib(k) and neib(L) don't intersect
    292296          if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.  ! a03 01
    293      #     -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then ! 06
    294 C Neighbours neib(k) and neib(L) cover fully the atom "i" and as(i)=0
     297     &     -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then ! 06
     298! Neighbours neib(k) and neib(L) cover fully the atom "i" and as(i)=0
    295299           As(i)=0.d0
    296300           goto 11
    297301          elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.  ! a03 02
    298      #      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
    299 C Don't exclude neib(k)
     302     &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
     303! Don't exclude neib(k)
    300304           neib(0)=neib(0)-1
    301305           do j=k,neib(0)
    … …  
    304308           goto 12
    305309          elseif(D.gt.0.d0) then                        ! a03 03
    306 C Circles neib(L) and neib(k) have two intersection points.
     310! Circles neib(L) and neib(k) have two intersection points.
    307311           am=2d0*((b1-b2)**2+(c1-c2)**2)
    308312           t=-2d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
    … …  
    324328
    325329         elseif(neibor(neib(k),1).lt.0d0.and.neibor(neib(L),1).gt.0d0)! a07
    326      #           then
     330     &           then
    327331          b1=neibor(neib(k),2)/neibor(neib(k),1)
    328332          c1=neibor(neib(k),3)/neibor(neib(k),1)
    … …  
    332336          d2=neibor(neib(L),4)/neibor(neib(L),1)
    333337          D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
    334      #        (d1-d2)**2)+(b1*c2-b2*c1)**2
    335 C If D<0 then the circles neib(k) and neib(L) don't intersect
     338     &        (d1-d2)**2)+(b1*c2-b2*c1)**2
     339! If D<0 then the circles neib(k) and neib(L) don't intersect
    336340          if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.   ! a07 01
    337      #     dsqrt(b2*b2+c2*c2-4d0*d2)-dsqrt(b1*b1+c1*c1-4d0*d1)) then    ! 10
    338 C atom "i" is covered fully by neib(k) and neib(L) and as(i)=0.
     341     &     dsqrt(b2*b2+c2*c2-4d0*d2)-dsqrt(b1*b1+c1*c1-4d0*d1)) then    ! 10
     342! atom "i" is covered fully by neib(k) and neib(L) and as(i)=0.
    339343           As(i)=0.d0
    340344           goto 12
    341345          elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.   ! a07 02
    342      #      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
    343 C discard the circle neib(L)
     346     &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
     347! discard the circle neib(L)
    344348           neib(0)=neib(0)-1
    345349           do j=L,neib(0)
    … …  
    347351           enddo
    348352          elseif(D.gt.0d0) then                          ! a07 03
    349 C There are two IP's between neib(k) and neib(L).
     353! There are two IP's between neib(k) and neib(L).
    350354           am=2d0*((b1-b2)**2+(c1-c2)**2)
    351355           t=-2d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
    … …  
    355359           p1=(c1-c2)*pom
    356360           p2=(b1-b2)*pom
    357 c-- Assign t and s coordinates and the order number of circles
     361!-- Assign t and s coordinates and the order number of circles
    358362           vertex(iv,1)=(t+p1)/am
    359363           vertex(iv,2)=(s-p2)/am
    … …  
    365369           vertex(iv,3)=neib(k)
    366370           vertex(iv,4)=neib(L)
    367 c-- 
     371!-- 
    368372          endif                                                 ! 10
    369373
    370374         elseif(neibor(neib(k),1).lt.0d0.and.neibor(neib(L),1).lt.0d0) ! a09
    371      #           then
     375     &           then
    372376          b1=neibor(neib(k),2)/neibor(neib(k),1)
    373377          c1=neibor(neib(k),3)/neibor(neib(k),1)
    … …  
    377381          d2=neibor(neib(L),4)/neibor(neib(L),1)
    378382          D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
    379      #        (d1-d2)**2)+(b1*c2-b2*c1)**2
    380 C D<0 - no intersection poit between neib(k) and neib(L)
     383     &        (d1-d2)**2)+(b1*c2-b2*c1)**2
     384! D<0 - no intersection poit between neib(k) and neib(L)
    381385          if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.
    382      #   dsqrt(b2*b2+c2*c2-4d0*d2)-dsqrt(b1*b1+c1*c1-4d0*d1)) then!12 ! a09 01
    383 C omit the circle neib(L)
     386     &   dsqrt(b2*b2+c2*c2-4d0*d2)-dsqrt(b1*b1+c1*c1-4d0*d1)) then!12 ! a09 01
     387! omit the circle neib(L)
    384388           neib(0)=neib(0)-1
    385389           do j=L,neib(0)
    … …  
    387391           enddo
    388392          elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le. ! a09 02
    389      #     -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
    390 C Omit the circle neib(k)
     393     &     -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
     394! Omit the circle neib(k)
    391395           neib(0)=neib(0)-1
    392396           do j=k,neib(0)
    … …  
    395399           goto 12
    396400          elseif(D.le.0.d0) then                            ! a09 03
    397 C The whole surface of atom "i" is covered fully, and as(i)=0.
     401! The whole surface of atom "i" is covered fully, and as(i)=0.
    398402           As(i)=0.d0
    399403           goto 11
    400404          else                                             ! a09 04
    401 C Two intersection points
     405! Two intersection points
    402406           am=2.d0*((b1-b2)**2+(c1-c2)**2)
    403407           t=-2.d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
    … …  
    407411           p1=(c1-c2)*pom
    408412           p2=(b1-b2)*pom
    409 C Assign the t and s coordinates and order numbers of the circles
     413! Assign the t and s coordinates and order numbers of the circles
    410414           vertex(iv,1)=(t+p1)/am
    411415           vertex(iv,2)=(s-p2)/am
    … …  
    428432        do j=1,neib(0)
    429433        if(neibor(neib(j),1).eq.0.d0) then
    430 C Collect the lines (after rotation it is empty).
     434! Collect the lines (after rotation it is empty).
    431435         ita2=ita2+1
    432436         ta2(ita2)=j
    433437        endif
    434438        if(neibor(neib(j),1).lt.0.d0) then
    435 C Collect the neighbours with inner part
     439! Collect the neighbours with inner part
    436440         ita3=ita3+1
    437441         ta3(ita3)=j
    … …  
    439443       enddo
    440444
    441 c*** Consider to remove this part.
     445!*** Consider to remove this part.
    442446         if(ita2.gt.0.and.ita3.lt.1) then
    443447           As(i)=R22p
    … …  
    472476          As(i)=0d0
    473477         endif
    474 c****** End of could-be-removed part
     478!****** End of could-be-removed part
    475479
    476480       neibp(0)=neib(0)
    … …  
    481485       enddo
    482486       do j=1,neib(0)
    483 C "iv" is the number of intersection point.
     487! "iv" is the number of intersection point.
    484488        do k=1,iv
    485489       if(neibp(j).eq.vertex(k,3).or.neibp(j).eq.vertex(k,4)) goto 30
    486490        enddo
    487491         ifu=ifu+1
    488 C Arrays with the ordering numbers of full arcs.
     492! Arrays with the ordering numbers of full arcs.
    489493         fullarc(ifu)=neibp(j)
    490494         al(ifu)=neibp(j)
    … …  
    500504       al(0)=ifu
    501505
    502 C Loop over full arcs.
     506! Loop over full arcs.
    503507       do k=1,ifu
    504508         a=neibor(fullarc(k),1)
    … …  
    512516         d=neibor(fullarc(k),4)
    513517         if(a.ne.0d0) then !True if rotation has taken place.
    514 C Remove the part of the atomic surface covered cut by this full arc
     518! Remove the part of the atomic surface covered cut by this full arc
    515519          As(i)=As(i)-aia*R22p*(1d0+(-d/a-R42)*dabs(a)/dsqrt((R42*a-
    516      #                 d)**2+R42*(b*b+c*c)))
     520     &                 d)**2+R42*(b*b+c*c)))
    517521          dadx(1,1)=2d0*ddat(fullarc(k),2)
    518522          dadx(1,2)=2d0*ddat(fullarc(k),3)
    … …  
    533537          gp(4)=(b*b+c*c-2d0*a*d+(R42+R42)*a*a)*div
    534538          grad(i,fullarc(k),1)=gp(1)*dadx(1,1)+gp(2)*dadx(2,1)+
    535      #        gp(3)*dadx(3,1)+gp(4)*dadx(4,1)
     539     &        gp(3)*dadx(3,1)+gp(4)*dadx(4,1)
    536540          grad(i,fullarc(k),2)=gp(1)*dadx(1,2)+gp(2)*dadx(2,2)+
    537      #        gp(3)*dadx(3,2)+gp(4)*dadx(4,2)
     541     &        gp(3)*dadx(3,2)+gp(4)*dadx(4,2)
    538542          grad(i,fullarc(k),3)=gp(1)*dadx(1,3)+gp(2)*dadx(2,3)+
    539      #        gp(3)*dadx(3,3)+gp(4)*dadx(4,3)
     543     &        gp(3)*dadx(3,3)+gp(4)*dadx(4,3)
    540544
    541545         else
    … …  
    544548       enddo
    545549
    546 c Loop over the circles with intersection points.
     550! Loop over the circles with intersection points.
    547551       do k=1,ine
    548552         jk=neib(k) ! The order number of kth neighbour.
    … …  
    602606             yv=(bcd+2d0*R42*ak)/daba*vv1
    603607             zv=wv/a*vv1
    604 C Move the (t,s) origo into the centre of this circle.
    605 C Modify t and s coordinates of IP.
     608! Move the (t,s) origo into the centre of this circle.
     609! Modify t and s coordinates of IP.
    606610          do j=1,nyx
    607611           ayx(j,1)=ayx(j,1)+ba ! t
    … …  
    612616          cs=-ca
    613617          rr=0.5d0*wv/daba
    614 C Calculate the polar angle of IP.
     618! Calculate the polar angle of IP.
    615619          do j=1,nyx
    616620           ayx1(j)=datan2(ayx(j,2),ayx(j,1))
    617621          enddo
    618 C Sort the angles in ascending order.
     622! Sort the angles in ascending order.
    619623          do j=1,nyx-1
    620624           a1=ayx1(j)
    … …  
    650654          do j=1,nyx
    651655
    652 c 
    653 c  Escape 'bad' intersections.
     656! 
     657!  Escape 'bad' intersections.
    654658           if(dabs(ayx1(j+1)-ayx1(j)).lt.1d-8) goto 40
    655659
    … …  
    657661           ap1=ct+rr*dcos(prob) ! The middle point of the arc.
    658662           ap2=cs+rr*dsin(prob)
    659 C Verify if the middle point belongs to a covered part.
    660 C If yes then omit.
     663! Verify if the middle point belongs to a covered part.
     664! If yes then omit.
    661665           do ij=1,ial
    662666          if(neibor(al(ij),1)*(ap1*ap1+ap2*ap2)+neibor(al(ij),2)*ap1+
    663      #         neibor(al(ij),3)*ap2+neibor(al(ij),4).lt.0d0) goto 40
     667     &         neibor(al(ij),3)*ap2+neibor(al(ij),4).lt.0d0) goto 40
    664668           enddo
    665669             alp=ayx1(j)
    … …  
    675679             dcotbma=1.d0/dtan(5d-1*(bet-alp))
    676680             uv=daba*(bcd+2d0*R42*ak)*dcbetmalp
    677      #          -a*dsqrt(b2c2-4d0*a*d)*(b*dcbetpalp+c*dsbetpalp)
     681     &          -a*dsqrt(b2c2-4d0*a*d)*(b*dcbetpalp+c*dsbetpalp)
    678682             As(i)=As(i)+R2*(aia*(alp-bet)+(d+R42*a)*vv1*(pi-2d0*
    679      #             datan(uv/(2d0*ak*vv*dsbetmalp))))
     683     &             datan(uv/(2d0*ak*vv*dsbetmalp))))
    680684             kk=vrx(j,4)
    681685             LL=vrx(j+1,4)
    … …  
    696700             am2=wv*aia*a2*(ab2*dsbet-ac2*dcbet)
    697701             dalp(1)=(b1*ab1+c1*ac1-2d0*d*a1*a1-a*
    698      1               (b1*b1+c1*c1-4d0*a1*d1)-2d0*d*aia*a1*(ab1*
    699      2               dcalp+ac1*dsalp)/wv+wv*aia*a1*
    700      3               (b1*dcalp+c1*dsalp))/am1
     702     &               (b1*b1+c1*c1-4d0*a1*d1)-2d0*d*aia*a1*(ab1*
     703     &               dcalp+ac1*dsalp)/wv+wv*aia*a1*
     704     &               (b1*dcalp+c1*dsalp))/am1
    701705             dalp(2)=(-a1*ab1+b*a1*a1+b*aia*a1*(ab1*
    702      1               dcalp+ac1*dsalp)/wv-wv*aia*a1*a1*dcalp)/am1
     706     &               dcalp+ac1*dsalp)/wv-wv*aia*a1*a1*dcalp)/am1
    703707             dalp(3)=(-a1*ac1+c*a1*a1+c*aia*a1*(ab1*
    704      #               dcalp+ac1*dsalp)/wv-wv*aia*a1*a1*dsalp)/am1
     708     &               dcalp+ac1*dsalp)/wv-wv*aia*a1*a1*dsalp)/am1
    705709             dalp(4)=(-2d0*a*a1*a1-2d0*daba*a1*(ab1*
    706      #               dcalp+ac1*dsalp)/wv)/am1
     710     &               dcalp+ac1*dsalp)/wv)/am1
    707711             dbet(1)=(b2*ab2+c2*ac2-2d0*d*a2*a2-a*
    708      #               (b2*b2+c2*c2-4d0*a2*d2)-2d0*d*aia*a2*(ab2*
    709      #               dcbet+ac2*dsbet)/wv+wv*aia*a2*
    710      #               (b2*dcbet+c2*dsbet))/am2
     712     &               (b2*b2+c2*c2-4d0*a2*d2)-2d0*d*aia*a2*(ab2*
     713     &               dcbet+ac2*dsbet)/wv+wv*aia*a2*
     714     &               (b2*dcbet+c2*dsbet))/am2
    711715             dbet(2)=(-a2*ab2+b*a2*a2+b*aia*a2*(ab2*
    712      #               dcbet+ac2*dsbet)/wv-wv*aia*a2*a2*dcbet)/am2
     716     &               dcbet+ac2*dsbet)/wv-wv*aia*a2*a2*dcbet)/am2
    713717             dbet(3)=(-a2*ac2+c*a2*a2+c*aia*a2*(ab2*
    714      #               dcbet+ac2*dsbet)/wv-wv*aia*a2*a2*dsbet)/am2
     718     &               dcbet+ac2*dsbet)/wv-wv*aia*a2*a2*dsbet)/am2
    715719             dbet(4)=(-2d0*a*a2*a2-2d0*daba*a2*(ab2*
    716      #               dcbet+ac2*dsbet)/wv)/am2
     720     &               dcbet+ac2*dsbet)/wv)/am2
    717721             daalp(1)=(-b*ab1-c*ac1+wv*wv*a1+2d0*ak*d1+
    718      #              wv*aia*((ab1-b*a1)*dcalp+(ac1-c*a1)*dsalp))/am1
     722     &              wv*aia*((ab1-b*a1)*dcalp+(ac1-c*a1)*dsalp))/am1
    719723             daalp(2)=(a*ab1-ak*b1+wv*daba*a1*dcalp)/am1
    720724             daalp(3)=(a*ac1-ak*c1+wv*daba*a1*dsalp)/am1
    721725             daalp(4)=(2d0*ak*a1)/am1
    722726             dabet(1)=(-b*ab2-c*ac2+wv*wv*a2+2d0*ak*d2+
    723      #             wv*aia*((ab2-b*a2)*dcbet+(ac2-c*a2)*dsbet))/am2
     727     &             wv*aia*((ab2-b*a2)*dcbet+(ac2-c*a2)*dsbet))/am2
    724728             dabet(2)=(a*ab2-ak*b2+wv*daba*a2*dcbet)/am2
    725729             dabet(3)=(a*ac2-ak*c2+wv*daba*a2*dsbet)/am2
    … …  
    734738             dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2     
    735739             dy(1)=(-2d0*d+4d0*R42*a)/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
    736      #             R42*ak)*(aia*vv+daba*dv(1))/ak*vv2
     740     &             R42*ak)*(aia*vv+daba*dv(1))/ak*vv2
    737741               dy(2)=2d0*b/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
    738      #             R42*ak)*aia*dv(2)/a*vv2
     742     &             R42*ak)*aia*dv(2)/a*vv2
    739743               dy(3)=2d0*c/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
    740      #             R42*ak)*aia*dv(3)/a*vv2
     744     &             R42*ak)*aia*dv(3)/a*vv2
    741745               dy(4)=-2.d0*a/daba*vv1-(b*b+c*c-2.d0*a*d+2.d0*
    742      #             R42*ak)*aia*dv(4)/a*vv2     
     746     &             R42*ak)*aia*dv(4)/a*vv2     
    743747             dz(1)=-2d0*d/wv/a*vv1-wv*(vv+a*dv(1))/ak*vv2
    744748             dz(2)=b/wv/a*vv1-wv*dv(2)/a*vv2
    … …  
    746750             dz(4)=-2d0/wv*vv1-wv*dv(4)/a*vv2
    747751             dt(1)=dcotbma*dy(1)-5d-1*yv*
    748      1               (dbet(1)-dalp(1))/dsbetmalp**2-((b*dcbetpalp+
    749      2             c*dsbetpalp)/dsbetmalp)*dz(1)-
    750      3             5d-1*zv*(((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp)*
    751      4             (dalp(1)+dbet(1))-((b*dcbetpalp+c*dsbetpalp)*
    752      5             dcbetmalp)*(dbet(1)-dalp(1)))/dsbetmalp**2
     752     &               (dbet(1)-dalp(1))/dsbetmalp**2-((b*dcbetpalp+
     753     &             c*dsbetpalp)/dsbetmalp)*dz(1)-
     754     &             5d-1*zv*(((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp)*
     755     &             (dalp(1)+dbet(1))-((b*dcbetpalp+c*dsbetpalp)*
     756     &             dcbetmalp)*(dbet(1)-dalp(1)))/dsbetmalp**2
    753757             dt(2)=dcotbma*dy(2)-5d-1*yv*
    754      1             (dbet(2)-dalp(2))/dsbetmalp**2-(b*dcbetpalp+
    755      2             c*dsbetpalp)/dsbetmalp*dz(2)-
    756      3             5d-1*zv*((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp*
    757      4             (dalp(2)+dbet(2))-(b*dcbetpalp+c*dsbetpalp)*
    758      5             dcbetmalp*(dbet(2)-dalp(2))+dcbetpalp
    759      6             *2d0*dsbetmalp)/dsbetmalp**2
     758     &             (dbet(2)-dalp(2))/dsbetmalp**2-(b*dcbetpalp+
     759     &             c*dsbetpalp)/dsbetmalp*dz(2)-
     760     &             5d-1*zv*((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp*
     761     &             (dalp(2)+dbet(2))-(b*dcbetpalp+c*dsbetpalp)*
     762     &             dcbetmalp*(dbet(2)-dalp(2))+dcbetpalp
     763     &             *2d0*dsbetmalp)/dsbetmalp**2
    760764             dt(3)=dcotbma*dy(3)-5d-1*yv*
    761      1             (dbet(3)-dalp(3))/dsbetmalp**2-(b*dcbetpalp+
    762      2             c*dsbetpalp)/dsbetmalp*dz(3)-
    763      3             5d-1*zv*((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp*
    764      4             (dalp(3)+dbet(3))-(b*dcbetpalp+c*dsbetpalp)*
    765      5             dcbetmalp*(dbet(3)-dalp(3))+dsbetpalp*2d0*
    766      6             dsbetmalp)/dsbetmalp**2
     765     &             (dbet(3)-dalp(3))/dsbetmalp**2-(b*dcbetpalp+
     766     &             c*dsbetpalp)/dsbetmalp*dz(3)-
     767     &             5d-1*zv*((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp*
     768     &             (dalp(3)+dbet(3))-(b*dcbetpalp+c*dsbetpalp)*
     769     &             dcbetmalp*(dbet(3)-dalp(3))+dsbetpalp*2d0*
     770     &             dsbetmalp)/dsbetmalp**2
    767771             dt(4)=dcotbma*dy(4)-5d-1*yv*
    768      1             (dbet(4)-dalp(4))/dsbetmalp**2-(b*dcbetpalp+
    769      2             c*dsbetpalp)/dsbetmalp*dz(4)-
    770      3             5d-1*zv*((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp*
    771      4             (dalp(4)+dbet(4))-(b*dcbetpalp+c*dsbetpalp)*
    772      5             dcbetmalp*(dbet(4)-dalp(4)))/dsbetmalp**2
     772     &             (dbet(4)-dalp(4))/dsbetmalp**2-(b*dcbetpalp+
     773     &             c*dsbetpalp)/dsbetmalp*dz(4)-
     774     &             5d-1*zv*((-b*dsbetpalp+c*dcbetpalp)*dsbetmalp*
     775     &             (dalp(4)+dbet(4))-(b*dcbetpalp+c*dsbetpalp)*
     776     &             dcbetmalp*(dbet(4)-dalp(4)))/dsbetmalp**2
    773777             di(1)=R2*(aia*(dalp(1)-dbet(1))+(pi-2d0*datan(5d-1*tt))*
    774      #               dx(1)-4d0*xv*dt(1)/(4d0+tt**2))
     778     &               dx(1)-4d0*xv*dt(1)/(4d0+tt**2))
    775779             di(2)=R2*(aia*(dalp(2)-dbet(2))+(pi-2d0*datan(5d-1*tt))*
    776      #               dx(2)-4d0*xv*dt(2)/(4d0+tt**2))
     780     &               dx(2)-4d0*xv*dt(2)/(4d0+tt**2))
    777781             di(3)=R2*(aia*(dalp(3)-dbet(3))+(pi-2d0*datan(5d-1*tt))*
    778      #               dx(3)-4d0*xv*dt(3)/(4d0+tt**2))
     782     &               dx(3)-4d0*xv*dt(3)/(4d0+tt**2))
    779783             di(4)=R2*(aia*(dalp(4)-dbet(4))+(pi-2d0*datan(5d-1*tt))*
    780      #               dx(4)-4d0*xv*dt(4)/(4d0+tt**2))
     784     &               dx(4)-4d0*xv*dt(4)/(4d0+tt**2))
    781785             dii(1,1)=2d0*ddat(jk,2)
    782786             dii(1,2)=2d0*ddat(jk,3)
    … …  
    794798             do idi=1,3
    795799              grad(i,jk,idi)=grad(i,jk,idi)+
    796      #                   di(1)*dii(1,idi)+di(2)*dii(2,idi)+
    797      #                   di(3)*dii(3,idi)+di(4)*dii(4,idi)
     800     &                   di(1)*dii(1,idi)+di(2)*dii(2,idi)+
     801     &                   di(3)*dii(3,idi)+di(4)*dii(4,idi)
    798802
    799803             enddo             
    800804             do idi=1,4
    801805               dta(idi)=5d-1*daalp(idi)*(((yv-zv*((-b*dsbetpalp
    802      1                  +c*dcbetpalp)*dsbetmalp+
    803      2                (b*dcbetpalp+c*dsbetpalp)*
    804      3                dcbetmalp))/dsbetmalp**2))
     806     &                  +c*dcbetpalp)*dsbetmalp+
     807     &                (b*dcbetpalp+c*dsbetpalp)*
     808     &                dcbetmalp))/dsbetmalp**2))
    805809               dtb(idi)=5d-1*dabet(idi)*(((-yv-zv*((-b*dsbetpalp
    806      1                  +c*dcbetpalp)*dsbetmalp-
    807      2                (b*dcbetpalp+c*dsbetpalp)*
    808      3                dcbetmalp))/dsbetmalp**2))
     810     &                  +c*dcbetpalp)*dsbetmalp-
     811     &                (b*dcbetpalp+c*dsbetpalp)*
     812     &                dcbetmalp))/dsbetmalp**2))
    809813             di1(idi)=R2*(aia*daalp(idi)-4d0*xv*dta(idi)/(4d0+tt*tt))
    810814            di2(idi)=R2*(-aia*dabet(idi)-4d0*xv*dtb(idi)/(4d0+tt*tt))
    … …  
    819823             do idi=1,3
    820824              grad(i,kk,idi)=grad(i,kk,idi)+
    821      1                   di1(1)*dii(1,idi)+di1(2)*dii(2,idi)+
    822      2                       di1(3)*dii(3,idi)+di1(4)*dii(4,idi)
     825     &                   di1(1)*dii(1,idi)+di1(2)*dii(2,idi)+
     826     &                       di1(3)*dii(3,idi)+di1(4)*dii(4,idi)
    823827             enddo
    824828
    … …  
    831835             do idi=1,3
    832836              grad(i,LL,idi)=grad(i,LL,idi)+
    833      1                   di2(1)*dii(1,idi)+di2(2)*dii(2,idi)+
    834      2                       di2(3)*dii(3,idi)+di2(4)*dii(4,idi)
     837     &                   di2(1)*dii(1,idi)+di2(2)*dii(2,idi)+
     838     &                       di2(3)*dii(3,idi)+di2(4)*dii(4,idi)
    835839             enddo
    836840
    … …  
    843847           a1=dsin(ang)
    844848           a2=dcos(ang)
    845 c Rotate the intersection points by the angle "ang".
    846 c After rotation the position of IP will be defined
    847 c by its first coordinate "ax(.,1)".
     849! Rotate the intersection points by the angle "ang".
     850! After rotation the position of IP will be defined
     851! by its first coordinate "ax(.,1)".
    848852           do j=1,nyx
    849853            ax(j,1)=ayx(j,1)*a2-ayx(j,2)*a1
    850854            ax(j,2)=ayx(j,1)*a1+ayx(j,2)*a2
    851855           enddo
    852 C Sorting by acsending order.
     856! Sorting by acsending order.
    853857           do j=1,nyx-1
    854858            a1=ax(j,1)
    … …  
    899903           p1=ayx(1,1)+c*probe(1)
    900904           p2=ayx(1,2)-b*probe(1)
    901 C Verify if the middle point belongs to a covered part.
    902 C Omit, if yes.
     905! Verify if the middle point belongs to a covered part.
     906! Omit, if yes.
    903907           do L=1,ial
    904908            if(neibor(al(L),1)*(p1*p1+p2*p2)+neibor(al(L),2)*p1+
    905      #         neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 20
     909     &         neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 20
    906910           enddo
    907911           arg2=bc*ayx1(1)+caba
    … …  
    914918            do L=1,ial
    915919             if(neibor(al(L),1)*(p1*p1+p2*p2)+neibor(al(L),2)*p1+
    916      #          neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 21
     920     &          neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 21
    917921            enddo
    918922            arg1=bc*ayx1(j-1)+caba
    … …  
    925929           do      L=1,ial
    926930            if(neibor(al(L),1)*(p1*p1+p2*p2)+neibor(al(L),2)*p1+
    927      1         neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 22
     931     &         neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 22
    928932           enddo
    929933           arg1=bc*ayx1(nyx)+caba
    … …  
    93493811      continue
    935939
    936 c
    937 c         changed in 2004
     940!
     941!         changed in 2004
    938942
    939943      do j=1,ine
    … …  
    952956
    953957      if(jjj.ne.0) then
    954 C Restore the original configuration if the molecule has been rotated.
    955 C This is necessary for calculation of gradients.
     958! Restore the original configuration if the molecule has been rotated.
     959! This is necessary for calculation of gradients.
    956960          xx=grad(i,i,1)
    957961          yy=grad(i,i,2)
    … …  
    988992111   continue
    989993
    990 c     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',     
    991 c    # '   gradz   Rad    Atom'
    992 c     write(3,*)
     994!     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',     
     995!    & '   gradz   Rad    Atom'
     996!     write(3,*)
    993997
    994998     
    … …  
    10001004        if (nmat(i)(1:1).eq.'h') ratom=0.0
    10011005
    1002 c     write(3,700),i,As(i),sigma(i),
    1003 c    #        enr,(gradan(i,k),k=1,3),ratom,nmat(i)
     1006!     write(3,700),i,As(i),sigma(i),
     1007!    &        enr,(gradan(i,k),k=1,3),ratom,nmat(i)
    10041008
    10051009      enddo
    10061010
    1007 c     write(3,*)'Total Area: ',sss
    1008 c     write(3,*)'Total solvation energy: ',energy
     1011!     write(3,*)'Total Area: ',sss
     1012!     write(3,*)'Total solvation energy: ',energy
    10091013
    10101014      esolan=sss
    10111015
    1012 c     close(3)
     1016!     close(3)
    10131017      return
    10141018
Note: See TracChangeset for help on using the changeset viewer.