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