| 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 | ! **************************************************************
|
|---|
| 11 |
|
|---|
| 12 | real*8 function esolan(nmol)
|
|---|
| 13 |
|
|---|
| 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
|
|---|
| 42 |
|
|---|
| 43 | include 'INCL.H'
|
|---|
| 44 |
|
|---|
| 45 |
|
|---|
| 46 | parameter (ks0=mxat*(mxat+1))
|
|---|
| 47 | parameter (ks2=mxat+mxat)
|
|---|
| 48 |
|
|---|
| 49 | dimension neib(0:mxat),vertex(ks0,4),ax(ks0,2),
|
|---|
| 50 | 1 pol(mxat),neibp(0:mxat),as(mxat),ayx(ks0,2),
|
|---|
| 51 | 2 ayx1(ks0),probe(ks0),dd(mxat),ddat(mxat,4),ivrx(ks0)
|
|---|
| 52 |
|
|---|
| 53 | dimension grad(mxat,mxat,3),dadx(4,3),gp(4),dalp(4),dbet(4),
|
|---|
| 54 | 1 daalp(4),dabet(4),vrx(ks0,4),dv(4),dx(4),dy(4),dz(4),dt(4),
|
|---|
| 55 | 2 di(4),dii(4,3),ss(mxat),dta(4),dtb(4),di1(4),di2(4),gs(3)
|
|---|
| 56 | dimension xold(-1:mxat),yold(-1:mxat),zold(-1:mxat)
|
|---|
| 57 | integer ta2(0:mxat),ta3(0:mxat),fullarc(0:mxat),al(0:ks2)
|
|---|
| 58 | real*8 neibor(mxat,4)
|
|---|
| 59 |
|
|---|
| 60 |
|
|---|
| 61 | ! open(3,file='solvation.dat')! Output file with solvation data
|
|---|
| 62 |
|
|---|
| 63 | energy=0.0d0 ! Total solvation energy
|
|---|
| 64 | sss=0.0d0 ! Total solvent accessible surface area
|
|---|
| 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
|
|---|
| 70 |
|
|---|
| 71 | do i=1,numat
|
|---|
| 72 | do j=1,3
|
|---|
| 73 | gradan(i,j)=0.0d0
|
|---|
| 74 | end do
|
|---|
| 75 | end do
|
|---|
| 76 |
|
|---|
| 77 | do i=1,numat
|
|---|
| 78 | xold(i)=xat(i)! Redefine the coordinates just for safety
|
|---|
| 79 | yold(i)=yat(i)!
|
|---|
| 80 | zold(i)=zat(i)!
|
|---|
| 81 | pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main'
|
|---|
| 82 | As(i)=pi4*pol(i)! Initially the whole surface of the atom is accessible.
|
|---|
| 83 | end do
|
|---|
| 84 |
|
|---|
| 85 | do iii=1,numat
|
|---|
| 86 | do jjj=1,numat
|
|---|
| 87 | grad(iii,jjj,1)=0.d0
|
|---|
| 88 | grad(iii,jjj,2)=0.d0
|
|---|
| 89 | grad(iii,jjj,3)=0.d0
|
|---|
| 90 | enddo
|
|---|
| 91 | enddo
|
|---|
| 92 |
|
|---|
| 93 | !***************************
|
|---|
| 94 | ! Start the computations *
|
|---|
| 95 | !***************************
|
|---|
| 96 |
|
|---|
| 97 | do 1 i=1,numat ! Lop over atoms "i"
|
|---|
| 98 | ! ----------------------------------------------------
|
|---|
| 99 |
|
|---|
| 100 | R=rvdw(i)
|
|---|
| 101 | jj=0
|
|---|
| 102 | do j=1,numat !Find the neighbours of "i"
|
|---|
| 103 | if(i.ne.j) then
|
|---|
| 104 | ddp=(xold(j)-xold(i))**2+(yold(j)-yold(i))**2
|
|---|
| 105 | & +(zold(j)-zold(i))**2! Distance between atom i and j
|
|---|
| 106 | if(ddp.lt.1e-10) then
|
|---|
| 107 | write(*,*)'ERROR in data: centres of two atoms coincide!'
|
|---|
| 108 | write(*,*)i,j,xold(i),yold(i),zold(i),rvdw(i),
|
|---|
| 109 | & xold(j),yold(j),zold(j),rvdw(j)
|
|---|
| 110 | stop 'Centres of atoms coincide!'
|
|---|
| 111 | endif
|
|---|
| 112 | ddp=dsqrt(ddp)
|
|---|
| 113 | if(ddp+R.le.rvdw(j)) then! Atom "i" is enclosed within "j"
|
|---|
| 114 | As(i)=0.d0 ! no accessible surface
|
|---|
| 115 | goto 1
|
|---|
| 116 | endif
|
|---|
| 117 | if(.not.(ddp.le.R-rvdw(j).or.ddp.ge.R+rvdw(j))) then
|
|---|
| 118 | jj=jj+1 ! The next neighbour of ith atom
|
|---|
| 119 | neib(jj)=j ! and its order number
|
|---|
| 120 | endif
|
|---|
| 121 | end if
|
|---|
| 122 | end do!Neighbours
|
|---|
| 123 |
|
|---|
| 124 | neib(0)=jj ! The number of neighbors of "i"
|
|---|
| 125 |
|
|---|
| 126 | if(neib(0).eq.0) goto 1 !Finish the atom i if it doesn't have neighbors
|
|---|
| 127 |
|
|---|
| 128 | !----
|
|---|
| 129 | R2=R*R ! R square
|
|---|
| 130 | R22=R2+R2 ! 2 * R square
|
|---|
| 131 | R42=R22+R22 ! 4 * R square
|
|---|
| 132 | R22p=R22*pi ! 2 pi R square
|
|---|
| 133 | R42p=R22p+R22p ! 4 pi R square
|
|---|
| 134 |
|
|---|
| 135 | do j=1,neib(0)
|
|---|
| 136 | k=neib(j)
|
|---|
| 137 | ddat(k,2)=xold(k)-xold(i)
|
|---|
| 138 | ddat(k,3)=yold(k)-yold(i)
|
|---|
| 139 | ddat(k,4)=zold(k)-zold(i)
|
|---|
| 140 | ddat(k,1)=rvdw(k)
|
|---|
| 141 | enddo
|
|---|
| 142 | ddat(i,2)=0d0
|
|---|
| 143 | ddat(i,3)=0d0
|
|---|
| 144 | ddat(i,4)=0d0
|
|---|
| 145 | ddat(i,1)=R
|
|---|
| 146 | !
|
|---|
| 147 | ! Verification of the north point of ith sphere
|
|---|
| 148 | !
|
|---|
| 149 | jjj=0
|
|---|
| 150 | !
|
|---|
| 151 | ! If jjj=0 then - no transformation
|
|---|
| 152 | ! else the transformation is necessary
|
|---|
| 153 | !
|
|---|
| 154 | k=neib(1) ! The order number of the first neighbour
|
|---|
| 155 | ! ddp - the square minimal distance of NP from the neighboring surfaces
|
|---|
| 156 | ddp=dabs((ddat(k,2))**2+(ddat(k,3))**2+
|
|---|
| 157 | & (R-ddat(k,4))**2-pol(k))
|
|---|
| 158 |
|
|---|
| 159 | do j=2,neib(0)
|
|---|
| 160 | k=neib(j)
|
|---|
| 161 | ddp2=dabs((ddat(k,2))**2+(ddat(k,3))**2+
|
|---|
| 162 | & (R-ddat(k,4))**2-pol(k))
|
|---|
| 163 | if(ddp2.lt.ddp) ddp=ddp2
|
|---|
| 164 | enddo
|
|---|
| 165 | !
|
|---|
| 166 | ! Check whether the NP of ith sphere is too close to the intersection line
|
|---|
| 167 | !
|
|---|
| 168 | do while (ddp.lt.1.d-6)
|
|---|
| 169 | jjj=jjj+1
|
|---|
| 170 | !
|
|---|
| 171 | ! generate grndom numbers
|
|---|
| 172 | !
|
|---|
| 173 | uga=grnd() ! Random \gamma angle
|
|---|
| 174 | ufi=grnd() ! Random \pi angle
|
|---|
| 175 | uga=pi*uga/2
|
|---|
| 176 | ufi=pi*ufi/2
|
|---|
| 177 | cf=dcos(ufi)
|
|---|
| 178 | sf=dsin(ufi)
|
|---|
| 179 | cg=dcos(uga)
|
|---|
| 180 | sg=dsin(uga)
|
|---|
| 181 | cfsg=cf*sg
|
|---|
| 182 | cfcg=cf*cg
|
|---|
| 183 | xx=R*cfsg
|
|---|
| 184 | yy=R*sf
|
|---|
| 185 | zz=R*cfcg
|
|---|
| 186 | k=neib(1)
|
|---|
| 187 | ddp=dabs((xx-ddat(k,2))**2+(yy-ddat(k,3))**2+
|
|---|
| 188 | & (zz-ddat(k,4))**2-pol(k))
|
|---|
| 189 | do j=2,neib(0)
|
|---|
| 190 | k=neib(j)
|
|---|
| 191 | ddp2=dabs((xx-ddat(k,2))**2+(yy-ddat(k,3))**2+
|
|---|
| 192 | & (zz-ddat(k,4))**2-pol(k))
|
|---|
| 193 | if(ddp2.lt.ddp) ddp=ddp2
|
|---|
| 194 | end do
|
|---|
| 195 | end do
|
|---|
| 196 | !
|
|---|
| 197 |
|
|---|
| 198 | if(jjj.ne.0) then ! Rotation is necessary
|
|---|
| 199 | sfsg=sf*sg
|
|---|
| 200 | sfcg=sf*cg
|
|---|
| 201 | do j=1,neib(0)
|
|---|
| 202 | k=neib(j)
|
|---|
| 203 | xx=ddat(k,2)
|
|---|
| 204 | yy=ddat(k,3)
|
|---|
| 205 | zz=ddat(k,4)
|
|---|
| 206 | ddat(k,2)=xx*cg-zz*sg ! (x) Coordinates
|
|---|
| 207 | ddat(k,3)=-xx*sfsg+yy*cf-zz*sfcg ! (y) after
|
|---|
| 208 | ddat(k,4)=xx*cfsg+yy*sf+zz*cfcg ! (z) rotation
|
|---|
| 209 | enddo
|
|---|
| 210 | endif
|
|---|
| 211 |
|
|---|
| 212 | ! iiiiiiiiiii
|
|---|
| 213 |
|
|---|
| 214 | pom=8.d0*pol(i)
|
|---|
| 215 |
|
|---|
| 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)
|
|---|
| 218 | do jj=1,neib(0)
|
|---|
| 219 | j=neib(jj)
|
|---|
| 220 | neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+
|
|---|
| 221 | & (ddat(j,4)-R)**2-pol(j) ! a
|
|---|
| 222 | neibor(j,2)=-pom*ddat(j,2) ! b
|
|---|
| 223 | neibor(j,3)=-pom*ddat(j,3) ! c
|
|---|
| 224 | neibor(j,4)=4d0*pol(i)*((ddat(j,2))**2+
|
|---|
| 225 | & (ddat(j,3))**2+(R+ddat(j,4))**2-pol(j)) ! d
|
|---|
| 226 | enddo
|
|---|
| 227 |
|
|---|
| 228 | ! end of the 1st cleaning
|
|---|
| 229 |
|
|---|
| 230 | iv=0
|
|---|
| 231 |
|
|---|
| 232 | nb=neib(0)
|
|---|
| 233 | k=neib(0)
|
|---|
| 234 | do while(k.gt.1) ! B
|
|---|
| 235 | k=k-1
|
|---|
| 236 | ! Analyse mutual disposition of every pair of neighbours
|
|---|
| 237 | do 13 L=neib(0),k+1,-1 ! A
|
|---|
| 238 |
|
|---|
| 239 | if(neibor(neib(k),1).gt.0d0.and.neibor(neib(L),1).gt.0d0) then ! 03 a01
|
|---|
| 240 |
|
|---|
| 241 | b1=neibor(neib(k),2)/neibor(neib(k),1)
|
|---|
| 242 | c1=neibor(neib(k),3)/neibor(neib(k),1)
|
|---|
| 243 | d1=neibor(neib(k),4)/neibor(neib(k),1)
|
|---|
| 244 | b2=neibor(neib(L),2)/neibor(neib(L),1)
|
|---|
| 245 | c2=neibor(neib(L),3)/neibor(neib(L),1)
|
|---|
| 246 | d2=neibor(neib(L),4)/neibor(neib(L),1)
|
|---|
| 247 | D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
|
|---|
| 248 | & (d1-d2)**2)+(b1*c2-b2*c1)**2
|
|---|
| 249 | ! if D<0 then the circles neib(k) and neib(L) don't intersect
|
|---|
| 250 | if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le. ! a01 01
|
|---|
| 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
|
|---|
| 253 | neib(0)=neib(0)-1
|
|---|
| 254 | do j=k,neib(0)
|
|---|
| 255 | neib(j)=neib(j+1)
|
|---|
| 256 | enddo
|
|---|
| 257 | goto 12
|
|---|
| 258 | elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le. ! a01 02
|
|---|
| 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
|
|---|
| 261 | neib(0)=neib(0)-1
|
|---|
| 262 | do j=L,neib(0)
|
|---|
| 263 | neib(j)=neib(j+1)
|
|---|
| 264 | enddo
|
|---|
| 265 | elseif(D.gt.0d0) then ! a01 03
|
|---|
| 266 | ! The circles nieb(L) and neib(k) have two intersection points (IP)
|
|---|
| 267 | am=2d0*((b1-b2)**2+(c1-c2)**2)
|
|---|
| 268 | t=-2d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
|
|---|
| 269 | s=-2d0*(c1-c2)*(d1-d2)+(b2*c1-b1*c2)*(b1-b2)
|
|---|
| 270 | iv=iv+1
|
|---|
| 271 | pom=dsqrt(D)
|
|---|
| 272 | p1=(c1-c2)*pom
|
|---|
| 273 | p2=(b1-b2)*pom
|
|---|
| 274 | vertex(iv,1)=(t+p1)/am ! t coordinate of the first IP
|
|---|
| 275 | vertex(iv,2)=(s-p2)/am ! s coordinate of the first IP
|
|---|
| 276 | vertex(iv,3)=neib(k) ! the order number of the first circle
|
|---|
| 277 | vertex(iv,4)=neib(L) ! the order number of the second circle
|
|---|
| 278 | iv=iv+1
|
|---|
| 279 | vertex(iv,1)=(t-p1)/am ! t coordinate of the second IP
|
|---|
| 280 | vertex(iv,2)=(s+p2)/am ! s coordinate of the second IP
|
|---|
| 281 | vertex(iv,3)=neib(k) ! the order number of the first circle
|
|---|
| 282 | vertex(iv,4)=neib(L) ! the order number of the second circle
|
|---|
| 283 | endif ! 04
|
|---|
| 284 |
|
|---|
| 285 | elseif(neibor(neib(k),1).gt.0d0.and. ! a03
|
|---|
| 286 | & neibor(neib(L),1).lt.0d0) then
|
|---|
| 287 | b1=neibor(neib(k),2)/neibor(neib(k),1)
|
|---|
| 288 | c1=neibor(neib(k),3)/neibor(neib(k),1)
|
|---|
| 289 | d1=neibor(neib(k),4)/neibor(neib(k),1)
|
|---|
| 290 | b2=neibor(neib(L),2)/neibor(neib(L),1)
|
|---|
| 291 | c2=neibor(neib(L),3)/neibor(neib(L),1)
|
|---|
| 292 | d2=neibor(neib(L),4)/neibor(neib(L),1)
|
|---|
| 293 | D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
|
|---|
| 294 | & (d1-d2)**2)+(b1*c2-b2*c1)**2
|
|---|
| 295 | ! if D<0 then neib(k) and neib(L) don't intersect
|
|---|
| 296 | if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le. ! a03 01
|
|---|
| 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
|
|---|
| 299 | As(i)=0.d0
|
|---|
| 300 | goto 11
|
|---|
| 301 | elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge. ! a03 02
|
|---|
| 302 | & dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
|
|---|
| 303 | ! Don't exclude neib(k)
|
|---|
| 304 | neib(0)=neib(0)-1
|
|---|
| 305 | do j=k,neib(0)
|
|---|
| 306 | neib(j)=neib(j+1)
|
|---|
| 307 | enddo
|
|---|
| 308 | goto 12
|
|---|
| 309 | elseif(D.gt.0.d0) then ! a03 03
|
|---|
| 310 | ! Circles neib(L) and neib(k) have two intersection points.
|
|---|
| 311 | am=2d0*((b1-b2)**2+(c1-c2)**2)
|
|---|
| 312 | t=-2d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
|
|---|
| 313 | s=-2d0*(c1-c2)*(d1-d2)+(b2*c1-b1*c2)*(b1-b2)
|
|---|
| 314 | iv=iv+1
|
|---|
| 315 | pom=dsqrt(D)
|
|---|
| 316 | p1=(c1-c2)*pom
|
|---|
| 317 | p2=(b1-b2)*pom
|
|---|
| 318 | vertex(iv,1)=(t+p1)/am ! t coordinate of the first IP
|
|---|
| 319 | vertex(iv,2)=(s-p2)/am ! s coordinate of the first IP
|
|---|
| 320 | vertex(iv,3)=neib(k) ! order number of the first circle
|
|---|
| 321 | vertex(iv,4)=neib(L) ! order number of the second circle
|
|---|
| 322 | iv=iv+1
|
|---|
| 323 | vertex(iv,1)=(t-p1)/am ! t coordinate of the second IP
|
|---|
| 324 | vertex(iv,2)=(s+p2)/am ! s coordinate of the second IP
|
|---|
| 325 | vertex(iv,3)=neib(k) ! order number of the first circle
|
|---|
| 326 | vertex(iv,4)=neib(L) ! order number of the second circle
|
|---|
| 327 | endif ! 06
|
|---|
| 328 |
|
|---|
| 329 | elseif(neibor(neib(k),1).lt.0d0.and.neibor(neib(L),1).gt.0d0)! a07
|
|---|
| 330 | & then
|
|---|
| 331 | b1=neibor(neib(k),2)/neibor(neib(k),1)
|
|---|
| 332 | c1=neibor(neib(k),3)/neibor(neib(k),1)
|
|---|
| 333 | d1=neibor(neib(k),4)/neibor(neib(k),1)
|
|---|
| 334 | b2=neibor(neib(L),2)/neibor(neib(L),1)
|
|---|
| 335 | c2=neibor(neib(L),3)/neibor(neib(L),1)
|
|---|
| 336 | d2=neibor(neib(L),4)/neibor(neib(L),1)
|
|---|
| 337 | D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
|
|---|
| 338 | & (d1-d2)**2)+(b1*c2-b2*c1)**2
|
|---|
| 339 | ! If D<0 then the circles neib(k) and neib(L) don't intersect
|
|---|
| 340 | if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le. ! a07 01
|
|---|
| 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.
|
|---|
| 343 | As(i)=0.d0
|
|---|
| 344 | goto 12
|
|---|
| 345 | elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge. ! a07 02
|
|---|
| 346 | & dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
|
|---|
| 347 | ! discard the circle neib(L)
|
|---|
| 348 | neib(0)=neib(0)-1
|
|---|
| 349 | do j=L,neib(0)
|
|---|
| 350 | neib(j)=neib(j+1)
|
|---|
| 351 | enddo
|
|---|
| 352 | elseif(D.gt.0d0) then ! a07 03
|
|---|
| 353 | ! There are two IP's between neib(k) and neib(L).
|
|---|
| 354 | am=2d0*((b1-b2)**2+(c1-c2)**2)
|
|---|
| 355 | t=-2d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
|
|---|
| 356 | s=-2d0*(c1-c2)*(d1-d2)+(b2*c1-b1*c2)*(b1-b2)
|
|---|
| 357 | iv=iv+1
|
|---|
| 358 | pom=dsqrt(D)
|
|---|
| 359 | p1=(c1-c2)*pom
|
|---|
| 360 | p2=(b1-b2)*pom
|
|---|
| 361 | !-- Assign t and s coordinates and the order number of circles
|
|---|
| 362 | vertex(iv,1)=(t+p1)/am
|
|---|
| 363 | vertex(iv,2)=(s-p2)/am
|
|---|
| 364 | vertex(iv,3)=neib(k)
|
|---|
| 365 | vertex(iv,4)=neib(L)
|
|---|
| 366 | iv=iv+1
|
|---|
| 367 | vertex(iv,1)=(t-p1)/am
|
|---|
| 368 | vertex(iv,2)=(s+p2)/am
|
|---|
| 369 | vertex(iv,3)=neib(k)
|
|---|
| 370 | vertex(iv,4)=neib(L)
|
|---|
| 371 | !--
|
|---|
| 372 | endif ! 10
|
|---|
| 373 |
|
|---|
| 374 | elseif(neibor(neib(k),1).lt.0d0.and.neibor(neib(L),1).lt.0d0) ! a09
|
|---|
| 375 | & then
|
|---|
| 376 | b1=neibor(neib(k),2)/neibor(neib(k),1)
|
|---|
| 377 | c1=neibor(neib(k),3)/neibor(neib(k),1)
|
|---|
| 378 | d1=neibor(neib(k),4)/neibor(neib(k),1)
|
|---|
| 379 | b2=neibor(neib(L),2)/neibor(neib(L),1)
|
|---|
| 380 | c2=neibor(neib(L),3)/neibor(neib(L),1)
|
|---|
| 381 | d2=neibor(neib(L),4)/neibor(neib(L),1)
|
|---|
| 382 | D=4d0*((b1-b2)*(b2*d1-b1*d2)+(c1-c2)*(c2*d1-c1*d2)-
|
|---|
| 383 | & (d1-d2)**2)+(b1*c2-b2*c1)**2
|
|---|
| 384 | ! D<0 - no intersection poit between neib(k) and neib(L)
|
|---|
| 385 | if(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.
|
|---|
| 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)
|
|---|
| 388 | neib(0)=neib(0)-1
|
|---|
| 389 | do j=L,neib(0)
|
|---|
| 390 | neib(j)=neib(j+1)
|
|---|
| 391 | enddo
|
|---|
| 392 | elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le. ! a09 02
|
|---|
| 393 | & -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
|
|---|
| 394 | ! Omit the circle neib(k)
|
|---|
| 395 | neib(0)=neib(0)-1
|
|---|
| 396 | do j=k,neib(0)
|
|---|
| 397 | neib(j)=neib(j+1)
|
|---|
| 398 | enddo
|
|---|
| 399 | goto 12
|
|---|
| 400 | elseif(D.le.0.d0) then ! a09 03
|
|---|
| 401 | ! The whole surface of atom "i" is covered fully, and as(i)=0.
|
|---|
| 402 | As(i)=0.d0
|
|---|
| 403 | goto 11
|
|---|
| 404 | else ! a09 04
|
|---|
| 405 | ! Two intersection points
|
|---|
| 406 | am=2.d0*((b1-b2)**2+(c1-c2)**2)
|
|---|
| 407 | t=-2.d0*(b1-b2)*(d1-d2)-(b2*c1-b1*c2)*(c1-c2)
|
|---|
| 408 | s=-2.d0*(c1-c2)*(d1-d2)+(b2*c1-b1*c2)*(b1-b2)
|
|---|
| 409 | iv=iv+1
|
|---|
| 410 | pom=dsqrt(D)
|
|---|
| 411 | p1=(c1-c2)*pom
|
|---|
| 412 | p2=(b1-b2)*pom
|
|---|
| 413 | ! Assign the t and s coordinates and order numbers of the circles
|
|---|
| 414 | vertex(iv,1)=(t+p1)/am
|
|---|
| 415 | vertex(iv,2)=(s-p2)/am
|
|---|
| 416 | vertex(iv,3)=neib(k)
|
|---|
| 417 | vertex(iv,4)=neib(L)
|
|---|
| 418 | iv=iv+1
|
|---|
| 419 | vertex(iv,1)=(t-p1)/am
|
|---|
| 420 | vertex(iv,2)=(s+p2)/am
|
|---|
| 421 | vertex(iv,3)=neib(k)
|
|---|
| 422 | vertex(iv,4)=neib(L)
|
|---|
| 423 | endif ! 12
|
|---|
| 424 |
|
|---|
| 425 | endif ! 03
|
|---|
| 426 | 13 continue
|
|---|
| 427 | 12 continue ! A
|
|---|
| 428 | enddo ! B
|
|---|
| 429 |
|
|---|
| 430 | ita2=0
|
|---|
| 431 | ita3=0
|
|---|
| 432 | do j=1,neib(0)
|
|---|
| 433 | if(neibor(neib(j),1).eq.0.d0) then
|
|---|
| 434 | ! Collect the lines (after rotation it is empty).
|
|---|
| 435 | ita2=ita2+1
|
|---|
| 436 | ta2(ita2)=j
|
|---|
| 437 | endif
|
|---|
| 438 | if(neibor(neib(j),1).lt.0.d0) then
|
|---|
| 439 | ! Collect the neighbours with inner part
|
|---|
| 440 | ita3=ita3+1
|
|---|
| 441 | ta3(ita3)=j
|
|---|
| 442 | endif
|
|---|
| 443 | enddo
|
|---|
| 444 |
|
|---|
| 445 | !*** Consider to remove this part.
|
|---|
| 446 | if(ita2.gt.0.and.ita3.lt.1) then
|
|---|
| 447 | As(i)=R22p
|
|---|
| 448 | if(ita2.gt.1) then
|
|---|
| 449 | amial=pi
|
|---|
| 450 | amaal=-pi
|
|---|
| 451 | amibe=pi
|
|---|
| 452 | amabe=-pi
|
|---|
| 453 | do j=1,ita2
|
|---|
| 454 | jj=ta2(j)
|
|---|
| 455 | p1=neibor(neib(jj),3)
|
|---|
| 456 | p2=neibor(neib(jj),2)
|
|---|
| 457 | alp=datan2(p1,p2)
|
|---|
| 458 | bet=alp
|
|---|
| 459 | if(alp.le.0d0) bet=bet+pi2
|
|---|
| 460 | if(amial.gt.alp) amial=alp
|
|---|
| 461 | if(amaal.lt.alp) amaal=alp
|
|---|
| 462 | if(amibe.gt.bet) amibe=bet
|
|---|
| 463 | if(amabe.lt.bet) amabe=bet
|
|---|
| 464 | enddo
|
|---|
| 465 | dal=amaal-amial
|
|---|
| 466 | dbe=amabe-amibe
|
|---|
| 467 | if(dal.lt.pi) then
|
|---|
| 468 | As(i)=R22*(pi-dal)
|
|---|
| 469 | elseif(dbe.lt.pi) then
|
|---|
| 470 | As(i)=R22*(pi-dbe)
|
|---|
| 471 | else
|
|---|
| 472 | As(i)=0d0
|
|---|
| 473 | endif
|
|---|
| 474 | endif
|
|---|
| 475 | elseif((nb.gt.0.and.neib(0).lt.1).or.ita3.gt.0) then
|
|---|
| 476 | As(i)=0d0
|
|---|
| 477 | endif
|
|---|
| 478 | !****** End of could-be-removed part
|
|---|
| 479 |
|
|---|
| 480 | neibp(0)=neib(0)
|
|---|
| 481 | ifu=0 ! Full arcs (without intersection points).
|
|---|
| 482 | ine=0 ! Circles with intersection points.
|
|---|
| 483 | do j=1,neib(0)
|
|---|
| 484 | neibp(j)=neib(j)
|
|---|
| 485 | enddo
|
|---|
| 486 | do j=1,neib(0)
|
|---|
| 487 | ! "iv" is the number of intersection point.
|
|---|
| 488 | do k=1,iv
|
|---|
| 489 | if(neibp(j).eq.vertex(k,3).or.neibp(j).eq.vertex(k,4)) goto 30
|
|---|
| 490 | enddo
|
|---|
| 491 | ifu=ifu+1
|
|---|
| 492 | ! Arrays with the ordering numbers of full arcs.
|
|---|
| 493 | fullarc(ifu)=neibp(j)
|
|---|
| 494 | al(ifu)=neibp(j)
|
|---|
| 495 | goto 31
|
|---|
| 496 | 30 continue
|
|---|
| 497 | ine=ine+1
|
|---|
| 498 | neib(ine)=neibp(j) !Array containing the intersection points.
|
|---|
| 499 | 31 continue
|
|---|
| 500 |
|
|---|
| 501 | enddo
|
|---|
| 502 | neib(0)=ine ! Number of IP.
|
|---|
| 503 | fullarc(0)=ifu ! Number of full arcs.
|
|---|
| 504 | al(0)=ifu
|
|---|
| 505 |
|
|---|
| 506 | ! Loop over full arcs.
|
|---|
| 507 | do k=1,ifu
|
|---|
| 508 | a=neibor(fullarc(k),1)
|
|---|
| 509 | if(a.lt.0.d0) then
|
|---|
| 510 | aia=-1.d0
|
|---|
| 511 | else
|
|---|
| 512 | aia=1.d0
|
|---|
| 513 | endif
|
|---|
| 514 | b=neibor(fullarc(k),2)
|
|---|
| 515 | c=neibor(fullarc(k),3)
|
|---|
| 516 | d=neibor(fullarc(k),4)
|
|---|
| 517 | if(a.ne.0d0) then !True if rotation has taken place.
|
|---|
| 518 | ! Remove the part of the atomic surface covered cut by this full arc
|
|---|
| 519 | As(i)=As(i)-aia*R22p*(1d0+(-d/a-R42)*dabs(a)/dsqrt((R42*a-
|
|---|
| 520 | & d)**2+R42*(b*b+c*c)))
|
|---|
| 521 | dadx(1,1)=2d0*ddat(fullarc(k),2)
|
|---|
| 522 | dadx(1,2)=2d0*ddat(fullarc(k),3)
|
|---|
| 523 | dadx(1,3)=2d0*(ddat(fullarc(k),4)-R)
|
|---|
| 524 | dadx(2,1)=-8d0*pol(i)
|
|---|
| 525 | dadx(2,2)=0d0
|
|---|
| 526 | dadx(2,3)=0d0
|
|---|
| 527 | dadx(3,1)=0d0
|
|---|
| 528 | dadx(3,2)=dadx(2,1)
|
|---|
| 529 | dadx(3,3)=0d0
|
|---|
| 530 | dadx(4,1)=8d0*pol(i)*ddat(fullarc(k),2)
|
|---|
| 531 | dadx(4,2)=8d0*pol(i)*ddat(fullarc(k),3)
|
|---|
| 532 | dadx(4,3)=8d0*pol(i)*(ddat(fullarc(k),4)+R)
|
|---|
| 533 | div=R22*R42p/((R42*a-d)**2+R42*(b*b+c*c))**(1.5)
|
|---|
| 534 | gp(1)=(R42*(b*b+c*c-2d0*a*d)+2d0*d*d)*div
|
|---|
| 535 | gp(2)=-b*(d+R42*a)*div
|
|---|
| 536 | gp(3)=-c*(d+R42*a)*div
|
|---|
| 537 | gp(4)=(b*b+c*c-2d0*a*d+(R42+R42)*a*a)*div
|
|---|
| 538 | grad(i,fullarc(k),1)=gp(1)*dadx(1,1)+gp(2)*dadx(2,1)+
|
|---|
| 539 | & gp(3)*dadx(3,1)+gp(4)*dadx(4,1)
|
|---|
| 540 | grad(i,fullarc(k),2)=gp(1)*dadx(1,2)+gp(2)*dadx(2,2)+
|
|---|
| 541 | & gp(3)*dadx(3,2)+gp(4)*dadx(4,2)
|
|---|
| 542 | grad(i,fullarc(k),3)=gp(1)*dadx(1,3)+gp(2)*dadx(2,3)+
|
|---|
| 543 | & gp(3)*dadx(3,3)+gp(4)*dadx(4,3)
|
|---|
| 544 |
|
|---|
| 545 | else
|
|---|
| 546 | As(i)=As(i)+R22p*d/dsqrt(R42*(b*b+c*c)+d*d)
|
|---|
| 547 | endif
|
|---|
| 548 | enddo
|
|---|
| 549 |
|
|---|
| 550 | ! Loop over the circles with intersection points.
|
|---|
| 551 | do k=1,ine
|
|---|
| 552 | jk=neib(k) ! The order number of kth neighbour.
|
|---|
| 553 | a=neibor(jk,1) ! a>0 - outer part, a<0-inner part
|
|---|
| 554 | ak=a*a
|
|---|
| 555 | daba=dabs(a)
|
|---|
| 556 | if(a.lt.0d0) then
|
|---|
| 557 | aia=-1d0
|
|---|
| 558 | else
|
|---|
| 559 | aia=1d0
|
|---|
| 560 | endif
|
|---|
| 561 | b=neibor(jk,2)
|
|---|
| 562 | c=neibor(jk,3)
|
|---|
| 563 | d=neibor(jk,4)
|
|---|
| 564 | b2c2=b*b+c*c
|
|---|
| 565 | vv=dsqrt((R42*a-d)**2+R42*b2c2)
|
|---|
| 566 | vv1=1d0/vv
|
|---|
| 567 | vv2=vv1*vv1
|
|---|
| 568 | wv=dsqrt(b2c2-4d0*a*d)
|
|---|
| 569 | nyx=0 ! The number of IP's which lie on this circle.
|
|---|
| 570 | ivr=0
|
|---|
| 571 | do j=1,iv
|
|---|
| 572 | if(vertex(j,3).eq.jk) then
|
|---|
| 573 | nyx=nyx+1
|
|---|
| 574 | ayx(nyx,1)=vertex(j,1) ! t coordinate of IP
|
|---|
| 575 | ayx(nyx,2)=vertex(j,2) ! s coordinate
|
|---|
| 576 | ivr=ivr+1
|
|---|
| 577 | ivrx(ivr)=j
|
|---|
| 578 | endif
|
|---|
| 579 | if(vertex(j,4).eq.jk) then
|
|---|
| 580 | nyx=nyx+1
|
|---|
| 581 | ayx(nyx,1)=vertex(j,1) ! t coordinate
|
|---|
| 582 | ayx(nyx,2)=vertex(j,2) ! s coordinate
|
|---|
| 583 | aa=vertex(j,3)
|
|---|
| 584 | vertex(j,3)=vertex(j,4)! make the fixed ordering number as the third
|
|---|
| 585 | vertex(j,4)=aa ! component in vertex if it was the 4th.
|
|---|
| 586 | ivr=ivr+1
|
|---|
| 587 | ivrx(ivr)=j
|
|---|
| 588 | endif
|
|---|
| 589 | enddo
|
|---|
| 590 |
|
|---|
| 591 | ial=ifu
|
|---|
| 592 | do j=1,ine
|
|---|
| 593 | if(neib(j).ne.jk) then
|
|---|
| 594 | ial=ial+1
|
|---|
| 595 | al(ial)=neib(j) ! Add other neighbours with IP's to this array
|
|---|
| 596 | endif
|
|---|
| 597 | enddo
|
|---|
| 598 | al(0)=ial ! The number of all "active" circles.
|
|---|
| 599 |
|
|---|
| 600 | if(a.ne.0.d0) then ! True after rotation
|
|---|
| 601 | aa=0.5d0/a
|
|---|
| 602 | ba=b*aa
|
|---|
| 603 | ca=c*aa
|
|---|
| 604 | bcd=b2c2-2d0*a*d
|
|---|
| 605 | xv=(d+R42*a)*vv1
|
|---|
| 606 | yv=(bcd+2d0*R42*ak)/daba*vv1
|
|---|
| 607 | zv=wv/a*vv1
|
|---|
| 608 | ! Move the (t,s) origo into the centre of this circle.
|
|---|
| 609 | ! Modify t and s coordinates of IP.
|
|---|
| 610 | do j=1,nyx
|
|---|
| 611 | ayx(j,1)=ayx(j,1)+ba ! t
|
|---|
| 612 | ayx(j,2)=ayx(j,2)+ca ! s
|
|---|
| 613 | enddo
|
|---|
| 614 |
|
|---|
| 615 | ct=-ba
|
|---|
| 616 | cs=-ca
|
|---|
| 617 | rr=0.5d0*wv/daba
|
|---|
| 618 | ! Calculate the polar angle of IP.
|
|---|
| 619 | do j=1,nyx
|
|---|
| 620 | ayx1(j)=datan2(ayx(j,2),ayx(j,1))
|
|---|
| 621 | enddo
|
|---|
| 622 | ! Sort the angles in ascending order.
|
|---|
| 623 | do j=1,nyx-1
|
|---|
| 624 | a1=ayx1(j)
|
|---|
| 625 | jj=j
|
|---|
| 626 | do ij=j+1,nyx
|
|---|
| 627 | if(a1.gt.ayx1(ij)) then
|
|---|
| 628 | a1=ayx1(ij)
|
|---|
| 629 | jj=ij
|
|---|
| 630 | endif
|
|---|
| 631 | enddo
|
|---|
| 632 | if(jj.ne.j) then
|
|---|
| 633 | a1=ayx1(jj)
|
|---|
| 634 | ayx1(jj)=ayx1(j)
|
|---|
| 635 | ayx1(j)=a1
|
|---|
| 636 | ivr=ivrx(jj)
|
|---|
| 637 | ivrx(jj)=ivrx(j)
|
|---|
| 638 | ivrx(j)=ivr
|
|---|
| 639 | endif
|
|---|
| 640 | enddo
|
|---|
| 641 |
|
|---|
| 642 | ayx1(nyx+1)=ayx1(1)+pi2
|
|---|
| 643 |
|
|---|
| 644 | do j=1,nyx
|
|---|
| 645 | do jkk=1,4
|
|---|
| 646 | vrx(j,jkk)=vertex(ivrx(j),jkk) ! vrx-helping array
|
|---|
| 647 | enddo
|
|---|
| 648 | enddo
|
|---|
| 649 |
|
|---|
| 650 | do jkk=1,4
|
|---|
| 651 | vrx(nyx+1,jkk)=vrx(1,jkk)
|
|---|
| 652 | enddo
|
|---|
| 653 |
|
|---|
| 654 | do j=1,nyx
|
|---|
| 655 |
|
|---|
| 656 | !
|
|---|
| 657 | ! Escape 'bad' intersections.
|
|---|
| 658 | if(dabs(ayx1(j+1)-ayx1(j)).lt.1d-8) goto 40
|
|---|
| 659 |
|
|---|
| 660 | prob=(ayx1(j)+ayx1(j+1))/2.d0
|
|---|
| 661 | ap1=ct+rr*dcos(prob) ! The middle point of the arc.
|
|---|
| 662 | ap2=cs+rr*dsin(prob)
|
|---|
| 663 | ! Verify if the middle point belongs to a covered part.
|
|---|
| 664 | ! If yes then omit.
|
|---|
| 665 | do ij=1,ial
|
|---|
| 666 | if(neibor(al(ij),1)*(ap1*ap1+ap2*ap2)+neibor(al(ij),2)*ap1+
|
|---|
| 667 | & neibor(al(ij),3)*ap2+neibor(al(ij),4).lt.0d0) goto 40
|
|---|
| 668 | enddo
|
|---|
| 669 | alp=ayx1(j)
|
|---|
| 670 | bet=ayx1(j+1)
|
|---|
| 671 | dsalp=dsin(alp)
|
|---|
| 672 | dcalp=dcos(alp)
|
|---|
| 673 | dsbet=dsin(bet)
|
|---|
| 674 | dcbet=dcos(bet)
|
|---|
| 675 | dcbetmalp=dcos(5d-1*(bet-alp))
|
|---|
| 676 | dcbetpalp=dcos(5d-1*(bet+alp))
|
|---|
| 677 | dsbetmalp=dsin(5d-1*(bet-alp))
|
|---|
| 678 | dsbetpalp=dsin(5d-1*(bet+alp))
|
|---|
| 679 | dcotbma=1.d0/dtan(5d-1*(bet-alp))
|
|---|
| 680 | uv=daba*(bcd+2d0*R42*ak)*dcbetmalp
|
|---|
| 681 | & -a*dsqrt(b2c2-4d0*a*d)*(b*dcbetpalp+c*dsbetpalp)
|
|---|
| 682 | As(i)=As(i)+R2*(aia*(alp-bet)+(d+R42*a)*vv1*(pi-2d0*
|
|---|
| 683 | & datan(uv/(2d0*ak*vv*dsbetmalp))))
|
|---|
| 684 | kk=vrx(j,4)
|
|---|
| 685 | LL=vrx(j+1,4)
|
|---|
| 686 | tt=yv*dcotbma-zv*(b*dcbetpalp+c*dsbetpalp)/dsbetmalp
|
|---|
| 687 | a1=neibor(kk,1)
|
|---|
| 688 | b1=neibor(kk,2)
|
|---|
| 689 | c1=neibor(kk,3)
|
|---|
| 690 | d1=neibor(kk,4)
|
|---|
| 691 | a2=neibor(LL,1)
|
|---|
| 692 | b2=neibor(LL,2)
|
|---|
| 693 | c2=neibor(LL,3)
|
|---|
| 694 | d2=neibor(LL,4)
|
|---|
| 695 | ab1=a*b1-b*a1
|
|---|
| 696 | ac1=a*c1-c*a1
|
|---|
| 697 | ab2=a*b2-b*a2
|
|---|
| 698 | ac2=a*c2-c*a2
|
|---|
| 699 | am1=wv*aia*a1*(ab1*dsalp-ac1*dcalp)
|
|---|
| 700 | am2=wv*aia*a2*(ab2*dsbet-ac2*dcbet)
|
|---|
| 701 | dalp(1)=(b1*ab1+c1*ac1-2d0*d*a1*a1-a*
|
|---|
| 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
|
|---|
| 705 | dalp(2)=(-a1*ab1+b*a1*a1+b*aia*a1*(ab1*
|
|---|
| 706 | & dcalp+ac1*dsalp)/wv-wv*aia*a1*a1*dcalp)/am1
|
|---|
| 707 | dalp(3)=(-a1*ac1+c*a1*a1+c*aia*a1*(ab1*
|
|---|
| 708 | & dcalp+ac1*dsalp)/wv-wv*aia*a1*a1*dsalp)/am1
|
|---|
| 709 | dalp(4)=(-2d0*a*a1*a1-2d0*daba*a1*(ab1*
|
|---|
| 710 | & dcalp+ac1*dsalp)/wv)/am1
|
|---|
| 711 | dbet(1)=(b2*ab2+c2*ac2-2d0*d*a2*a2-a*
|
|---|
| 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
|
|---|
| 715 | dbet(2)=(-a2*ab2+b*a2*a2+b*aia*a2*(ab2*
|
|---|
| 716 | & dcbet+ac2*dsbet)/wv-wv*aia*a2*a2*dcbet)/am2
|
|---|
| 717 | dbet(3)=(-a2*ac2+c*a2*a2+c*aia*a2*(ab2*
|
|---|
| 718 | & dcbet+ac2*dsbet)/wv-wv*aia*a2*a2*dsbet)/am2
|
|---|
| 719 | dbet(4)=(-2d0*a*a2*a2-2d0*daba*a2*(ab2*
|
|---|
| 720 | & dcbet+ac2*dsbet)/wv)/am2
|
|---|
| 721 | daalp(1)=(-b*ab1-c*ac1+wv*wv*a1+2d0*ak*d1+
|
|---|
| 722 | & wv*aia*((ab1-b*a1)*dcalp+(ac1-c*a1)*dsalp))/am1
|
|---|
| 723 | daalp(2)=(a*ab1-ak*b1+wv*daba*a1*dcalp)/am1
|
|---|
| 724 | daalp(3)=(a*ac1-ak*c1+wv*daba*a1*dsalp)/am1
|
|---|
| 725 | daalp(4)=(2d0*ak*a1)/am1
|
|---|
| 726 | dabet(1)=(-b*ab2-c*ac2+wv*wv*a2+2d0*ak*d2+
|
|---|
| 727 | & wv*aia*((ab2-b*a2)*dcbet+(ac2-c*a2)*dsbet))/am2
|
|---|
| 728 | dabet(2)=(a*ab2-ak*b2+wv*daba*a2*dcbet)/am2
|
|---|
| 729 | dabet(3)=(a*ac2-ak*c2+wv*daba*a2*dsbet)/am2
|
|---|
| 730 | dabet(4)=(2d0*ak*a2)/am2
|
|---|
| 731 | dv(2)=R42*b*vv1
|
|---|
| 732 | dv(3)=R42*c*vv1
|
|---|
| 733 | dv(4)=(d-R42*a)*vv1
|
|---|
| 734 | dv(1)=-dv(4)*R42
|
|---|
| 735 | dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2
|
|---|
| 736 | dx(2)=-(d+R42*a)*dv(2)*vv2
|
|---|
| 737 | dx(3)=-(d+R42*a)*dv(3)*vv2
|
|---|
| 738 | dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2
|
|---|
| 739 | dy(1)=(-2d0*d+4d0*R42*a)/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
|
|---|
| 740 | & R42*ak)*(aia*vv+daba*dv(1))/ak*vv2
|
|---|
| 741 | dy(2)=2d0*b/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
|
|---|
| 742 | & R42*ak)*aia*dv(2)/a*vv2
|
|---|
| 743 | dy(3)=2d0*c/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
|
|---|
| 744 | & R42*ak)*aia*dv(3)/a*vv2
|
|---|
| 745 | dy(4)=-2.d0*a/daba*vv1-(b*b+c*c-2.d0*a*d+2.d0*
|
|---|
| 746 | & R42*ak)*aia*dv(4)/a*vv2
|
|---|
| 747 | dz(1)=-2d0*d/wv/a*vv1-wv*(vv+a*dv(1))/ak*vv2
|
|---|
| 748 | dz(2)=b/wv/a*vv1-wv*dv(2)/a*vv2
|
|---|
| 749 | dz(3)=c/wv/a*vv1-wv*dv(3)/a*vv2
|
|---|
| 750 | dz(4)=-2d0/wv*vv1-wv*dv(4)/a*vv2
|
|---|
| 751 | dt(1)=dcotbma*dy(1)-5d-1*yv*
|
|---|
| 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
|
|---|
| 757 | dt(2)=dcotbma*dy(2)-5d-1*yv*
|
|---|
| 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
|
|---|
| 764 | dt(3)=dcotbma*dy(3)-5d-1*yv*
|
|---|
| 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
|
|---|
| 771 | dt(4)=dcotbma*dy(4)-5d-1*yv*
|
|---|
| 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
|
|---|
| 777 | di(1)=R2*(aia*(dalp(1)-dbet(1))+(pi-2d0*datan(5d-1*tt))*
|
|---|
| 778 | & dx(1)-4d0*xv*dt(1)/(4d0+tt**2))
|
|---|
| 779 | di(2)=R2*(aia*(dalp(2)-dbet(2))+(pi-2d0*datan(5d-1*tt))*
|
|---|
| 780 | & dx(2)-4d0*xv*dt(2)/(4d0+tt**2))
|
|---|
| 781 | di(3)=R2*(aia*(dalp(3)-dbet(3))+(pi-2d0*datan(5d-1*tt))*
|
|---|
| 782 | & dx(3)-4d0*xv*dt(3)/(4d0+tt**2))
|
|---|
| 783 | di(4)=R2*(aia*(dalp(4)-dbet(4))+(pi-2d0*datan(5d-1*tt))*
|
|---|
| 784 | & dx(4)-4d0*xv*dt(4)/(4d0+tt**2))
|
|---|
| 785 | dii(1,1)=2d0*ddat(jk,2)
|
|---|
| 786 | dii(1,2)=2d0*ddat(jk,3)
|
|---|
| 787 | dii(1,3)=2d0*(ddat(jk,4)-R)
|
|---|
| 788 | dii(2,1)=-R42-R42
|
|---|
| 789 | dii(2,2)=0d0
|
|---|
| 790 | dii(2,3)=0d0
|
|---|
| 791 | dii(3,1)=0d0
|
|---|
| 792 | dii(3,2)=dii(2,1)
|
|---|
| 793 | dii(3,3)=0d0
|
|---|
| 794 | dii(4,1)=dii(1,1)*R42
|
|---|
| 795 | dii(4,2)=dii(1,2)*R42
|
|---|
| 796 | dii(4,3)=2d0*(ddat(jk,4)+R)*R42
|
|---|
| 797 |
|
|---|
| 798 | do idi=1,3
|
|---|
| 799 | grad(i,jk,idi)=grad(i,jk,idi)+
|
|---|
| 800 | & di(1)*dii(1,idi)+di(2)*dii(2,idi)+
|
|---|
| 801 | & di(3)*dii(3,idi)+di(4)*dii(4,idi)
|
|---|
| 802 |
|
|---|
| 803 | enddo
|
|---|
| 804 | do idi=1,4
|
|---|
| 805 | dta(idi)=5d-1*daalp(idi)*(((yv-zv*((-b*dsbetpalp
|
|---|
| 806 | & +c*dcbetpalp)*dsbetmalp+
|
|---|
| 807 | & (b*dcbetpalp+c*dsbetpalp)*
|
|---|
| 808 | & dcbetmalp))/dsbetmalp**2))
|
|---|
| 809 | dtb(idi)=5d-1*dabet(idi)*(((-yv-zv*((-b*dsbetpalp
|
|---|
| 810 | & +c*dcbetpalp)*dsbetmalp-
|
|---|
| 811 | & (b*dcbetpalp+c*dsbetpalp)*
|
|---|
| 812 | & dcbetmalp))/dsbetmalp**2))
|
|---|
| 813 | di1(idi)=R2*(aia*daalp(idi)-4d0*xv*dta(idi)/(4d0+tt*tt))
|
|---|
| 814 | di2(idi)=R2*(-aia*dabet(idi)-4d0*xv*dtb(idi)/(4d0+tt*tt))
|
|---|
| 815 | enddo
|
|---|
| 816 |
|
|---|
| 817 | dii(1,1)=2d0*ddat(kk,2)
|
|---|
| 818 | dii(1,2)=2d0*ddat(kk,3)
|
|---|
| 819 | dii(1,3)=2d0*(ddat(kk,4)-R)
|
|---|
| 820 | dii(4,1)=dii(1,1)*R42
|
|---|
| 821 | dii(4,2)=dii(1,2)*R42
|
|---|
| 822 | dii(4,3)=2d0*(ddat(kk,4)+R)*R42
|
|---|
| 823 | do idi=1,3
|
|---|
| 824 | grad(i,kk,idi)=grad(i,kk,idi)+
|
|---|
| 825 | & di1(1)*dii(1,idi)+di1(2)*dii(2,idi)+
|
|---|
| 826 | & di1(3)*dii(3,idi)+di1(4)*dii(4,idi)
|
|---|
| 827 | enddo
|
|---|
| 828 |
|
|---|
| 829 | dii(1,1)=2d0*ddat(LL,2)
|
|---|
| 830 | dii(1,2)=2d0*ddat(LL,3)
|
|---|
| 831 | dii(1,3)=2d0*(ddat(LL,4)-R)
|
|---|
| 832 | dii(4,1)=dii(1,1)*R42
|
|---|
| 833 | dii(4,2)=dii(1,2)*R42
|
|---|
| 834 | dii(4,3)=2d0*(ddat(LL,4)+R)*R42
|
|---|
| 835 | do idi=1,3
|
|---|
| 836 | grad(i,LL,idi)=grad(i,LL,idi)+
|
|---|
| 837 | & di2(1)*dii(1,idi)+di2(2)*dii(2,idi)+
|
|---|
| 838 | & di2(3)*dii(3,idi)+di2(4)*dii(4,idi)
|
|---|
| 839 | enddo
|
|---|
| 840 |
|
|---|
| 841 | 40 continue
|
|---|
| 842 | enddo
|
|---|
| 843 |
|
|---|
| 844 | else ! analyse the case with lines (not necessary after rotation).
|
|---|
| 845 |
|
|---|
| 846 | ang=datan2(b,c)
|
|---|
| 847 | a1=dsin(ang)
|
|---|
| 848 | a2=dcos(ang)
|
|---|
| 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)".
|
|---|
| 852 | do j=1,nyx
|
|---|
| 853 | ax(j,1)=ayx(j,1)*a2-ayx(j,2)*a1
|
|---|
| 854 | ax(j,2)=ayx(j,1)*a1+ayx(j,2)*a2
|
|---|
| 855 | enddo
|
|---|
| 856 | ! Sorting by acsending order.
|
|---|
| 857 | do j=1,nyx-1
|
|---|
| 858 | a1=ax(j,1)
|
|---|
| 859 | jj=j
|
|---|
| 860 | do L=j+1,nyx
|
|---|
| 861 | if(a1.gt.ax(L,1)) then
|
|---|
| 862 | a1=ax(L,1)
|
|---|
| 863 | jj=L
|
|---|
| 864 | endif
|
|---|
| 865 | enddo
|
|---|
| 866 | if(jj.ne.j) then
|
|---|
| 867 | a1=ax(jj,1)
|
|---|
| 868 | a2=ax(jj,2)
|
|---|
| 869 | ax(jj,1)=ax(j,1)
|
|---|
| 870 | ax(jj,2)=ax(j,2)
|
|---|
| 871 | ax(j,1)=a1
|
|---|
| 872 | ax(j,2)=a2
|
|---|
| 873 | a1=ayx(jj,1)
|
|---|
| 874 | a2=ayx(jj,2)
|
|---|
| 875 | ayx(jj,1)=ayx(j,1)
|
|---|
| 876 | ayx(jj,2)=ayx(j,2)
|
|---|
| 877 | ayx(j,1)=a1
|
|---|
| 878 | ayx(j,2)=a2
|
|---|
| 879 | endif
|
|---|
| 880 | enddo
|
|---|
| 881 |
|
|---|
| 882 | if(c.ne.0d0) then
|
|---|
| 883 | do j=1,nyx
|
|---|
| 884 | ayx1(j)=(ayx(j,1)-ayx(1,1))/c
|
|---|
| 885 | enddo
|
|---|
| 886 | else
|
|---|
| 887 | do j=1,nyx
|
|---|
| 888 | ayx1(j)=-(ayx(j,2)-ayx(1,2))/b
|
|---|
| 889 | enddo
|
|---|
| 890 | endif
|
|---|
| 891 |
|
|---|
| 892 | probe(1)=-1d0
|
|---|
| 893 | do j=1,nyx-1
|
|---|
| 894 | probe(j+1)=(ayx1(j)+ayx1(j+1))/2d0
|
|---|
| 895 | enddo
|
|---|
| 896 | probe(nyx+1)=ayx1(nyx)+1d0
|
|---|
| 897 |
|
|---|
| 898 |
|
|---|
| 899 | bc=b*b+c*c
|
|---|
| 900 | cc=dsqrt(R42*(bc)+d*d)
|
|---|
| 901 | caba=c*ayx(1,1)-b*ayx(1,2)
|
|---|
| 902 |
|
|---|
| 903 | p1=ayx(1,1)+c*probe(1)
|
|---|
| 904 | p2=ayx(1,2)-b*probe(1)
|
|---|
| 905 | ! Verify if the middle point belongs to a covered part.
|
|---|
| 906 | ! Omit, if yes.
|
|---|
| 907 | do L=1,ial
|
|---|
| 908 | if(neibor(al(L),1)*(p1*p1+p2*p2)+neibor(al(L),2)*p1+
|
|---|
| 909 | & neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 20
|
|---|
| 910 | enddo
|
|---|
| 911 | arg2=bc*ayx1(1)+caba
|
|---|
| 912 | As(i)=As(i)+R22*d*(datan(arg2/cc)+pi/2d0)/cc
|
|---|
| 913 | 20 continue
|
|---|
| 914 |
|
|---|
| 915 | do j=2,nyx
|
|---|
| 916 | p1=ayx(1,1)+c*probe(j)
|
|---|
| 917 | p2=ayx(1,2)-b*probe(j)
|
|---|
| 918 | do L=1,ial
|
|---|
| 919 | if(neibor(al(L),1)*(p1*p1+p2*p2)+neibor(al(L),2)*p1+
|
|---|
| 920 | & neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 21
|
|---|
| 921 | enddo
|
|---|
| 922 | arg1=bc*ayx1(j-1)+caba
|
|---|
| 923 | arg2=bc*ayx1(j)+caba
|
|---|
| 924 | As(i)=As(i)+R22*d*(datan(arg2/cc)-datan(arg1/cc))/cc
|
|---|
| 925 | 21 continue
|
|---|
| 926 | enddo
|
|---|
| 927 | p1=ayx(1,1)+c*probe(nyx+1)
|
|---|
| 928 | p2=ayx(1,2)-b*probe(nyx+1)
|
|---|
| 929 | do L=1,ial
|
|---|
| 930 | if(neibor(al(L),1)*(p1*p1+p2*p2)+neibor(al(L),2)*p1+
|
|---|
| 931 | & neibor(al(L),3)*p2+neibor(al(L),4).lt.0d0) goto 22
|
|---|
| 932 | enddo
|
|---|
| 933 | arg1=bc*ayx1(nyx)+caba
|
|---|
| 934 | As(i)=As(i)+R22*d*(pi/2d0-datan(arg1/cc))/cc
|
|---|
| 935 | 22 continue
|
|---|
| 936 | endif
|
|---|
| 937 | enddo
|
|---|
| 938 | 11 continue
|
|---|
| 939 |
|
|---|
| 940 | !
|
|---|
| 941 | ! changed in 2004
|
|---|
| 942 |
|
|---|
| 943 | do j=1,ine
|
|---|
| 944 | al(ifu+j)=neib(j)
|
|---|
| 945 | enddo
|
|---|
| 946 | al(0)=ifu+ine
|
|---|
| 947 |
|
|---|
| 948 | do k=1,3
|
|---|
| 949 | ss(k)=0d0
|
|---|
| 950 | do j2=1,al(0)
|
|---|
| 951 | j=al(j2)
|
|---|
| 952 | ss(k)=ss(k)+grad(i,j,k)
|
|---|
| 953 | enddo
|
|---|
| 954 | grad(i,i,k)=-ss(k)
|
|---|
| 955 | enddo
|
|---|
| 956 |
|
|---|
| 957 | if(jjj.ne.0) then
|
|---|
| 958 | ! Restore the original configuration if the molecule has been rotated.
|
|---|
| 959 | ! This is necessary for calculation of gradients.
|
|---|
| 960 | xx=grad(i,i,1)
|
|---|
| 961 | yy=grad(i,i,2)
|
|---|
| 962 | zz=grad(i,i,3)
|
|---|
| 963 | grad(i,i,1)=xx*cg-yy*sfsg+zz*cfsg
|
|---|
| 964 | grad(i,i,2)=yy*cf+zz*sf
|
|---|
| 965 | grad(i,i,3)=-xx*sg-yy*sfcg+zz*cfcg
|
|---|
| 966 |
|
|---|
| 967 | do j2=1,al(0)
|
|---|
| 968 | j=al(j2)
|
|---|
| 969 | xx=grad(i,j,1)
|
|---|
| 970 | yy=grad(i,j,2)
|
|---|
| 971 | zz=grad(i,j,3)
|
|---|
| 972 | grad(i,j,1)=xx*cg-yy*sfsg+zz*cfsg
|
|---|
| 973 | grad(i,j,2)=yy*cf+zz*sf
|
|---|
| 974 | grad(i,j,3)=-xx*sg-yy*sfcg+zz*cfcg
|
|---|
| 975 | enddo
|
|---|
| 976 |
|
|---|
| 977 | endif
|
|---|
| 978 |
|
|---|
| 979 |
|
|---|
| 980 | 1 sss=sss+As(i)*sigma(i)
|
|---|
| 981 |
|
|---|
| 982 | do j=1,numat
|
|---|
| 983 | do k=1,3
|
|---|
| 984 | gs(k)=0.d0
|
|---|
| 985 | do i=1,numat
|
|---|
| 986 | gs(k)=gs(k)+sigma(i)*grad(i,j,k)
|
|---|
| 987 | enddo
|
|---|
| 988 | gradan(j,k)=gs(k)
|
|---|
| 989 | enddo
|
|---|
| 990 | enddo
|
|---|
| 991 |
|
|---|
| 992 | 111 continue
|
|---|
| 993 |
|
|---|
| 994 | ! write(3,*)' No Area sigma Enrg gradx grady',
|
|---|
| 995 | ! & ' gradz Rad Atom'
|
|---|
| 996 | ! write(3,*)
|
|---|
| 997 |
|
|---|
| 998 |
|
|---|
| 999 | do i=1,numat
|
|---|
| 1000 | enr=As(i)*sigma(i)
|
|---|
| 1001 | energy=energy+enr
|
|---|
| 1002 | ratom=rvdw(i)-rwater
|
|---|
| 1003 |
|
|---|
| 1004 | if (nmat(i)(1:1).eq.'h') ratom=0.0
|
|---|
| 1005 |
|
|---|
| 1006 | ! write(3,700),i,As(i),sigma(i),
|
|---|
| 1007 | ! & enr,(gradan(i,k),k=1,3),ratom,nmat(i)
|
|---|
| 1008 |
|
|---|
| 1009 | enddo
|
|---|
| 1010 |
|
|---|
| 1011 | ! write(3,*)'Total Area: ',sss
|
|---|
| 1012 | ! write(3,*)'Total solvation energy: ',energy
|
|---|
| 1013 |
|
|---|
| 1014 | esolan=sss
|
|---|
| 1015 |
|
|---|
| 1016 | ! close(3)
|
|---|
| 1017 | return
|
|---|
| 1018 |
|
|---|
| 1019 | 100 format(20i4)
|
|---|
| 1020 | 200 format(2i4,3f16.6)
|
|---|
| 1021 | 700 format(i5,7f8.3,5x,a4)
|
|---|
| 1022 |
|
|---|
| 1023 | end
|
|---|
