1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  2003  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     #################################################################
9c     ##                                                             ##
10c     ##  subroutine epitors2  --  pi-system torsion Hessian; numer  ##
11c     ##                                                             ##
12c     #################################################################
13c
14c
15c     "epitors2" calculates the second derivatives of the pi-system
16c     torsion energy for a single atom using finite difference methods
17c
18c
19      subroutine epitors2 (i)
20      use angbnd
21      use atoms
22      use group
23      use hessn
24      use pitors
25      use usage
26      implicit none
27      integer i,j,ipitors
28      integer ia,ib,ic
29      integer id,ie,ig
30      real*8 eps,fgrp
31      real*8 old,term
32      real*8, allocatable :: de(:,:)
33      real*8, allocatable :: d0(:,:)
34      logical proceed
35      logical twosided
36c
37c
38c     set stepsize for derivatives and default group weight
39c
40      eps = 1.0d-5
41      fgrp = 1.0d0
42      twosided = .false.
43      if (n .le. 50)  twosided = .true.
44c
45c     perform dynamic allocation of some local arrays
46c
47      allocate (de(3,n))
48      allocate (d0(3,n))
49c
50c     calculate numerical pi-system torsion Hessian for current atom
51c
52      do ipitors = 1, npitors
53         ia = ipit(1,ipitors)
54         ib = ipit(2,ipitors)
55         ic = ipit(3,ipitors)
56         id = ipit(4,ipitors)
57         ie = ipit(5,ipitors)
58         ig = ipit(6,ipitors)
59c
60c     decide whether to compute the current interaction
61c
62         proceed = .true.
63         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,id,ie,ig)
64         if (proceed)  proceed = (use(ia) .or. use(ib) .or. use(ic) .or.
65     &                              use(id) .or. use(ie) .or. use(ig))
66         if (proceed) then
67            term = fgrp / eps
68c
69c     find first derivatives for the base structure
70c
71            if (.not. twosided) then
72               call epitors2a (ipitors,de)
73               do j = 1, 3
74                  d0(j,ia) = de(j,ia)
75                  d0(j,ib) = de(j,ib)
76                  d0(j,ic) = de(j,ic)
77                  d0(j,id) = de(j,id)
78                  d0(j,ie) = de(j,ie)
79                  d0(j,ig) = de(j,ig)
80               end do
81            end if
82c
83c     find numerical x-components via perturbed structures
84c
85            old = x(i)
86            if (twosided) then
87               x(i) = x(i) - 0.5d0*eps
88               call epitors2a (ipitors,de)
89               do j = 1, 3
90                  d0(j,ia) = de(j,ia)
91                  d0(j,ib) = de(j,ib)
92                  d0(j,ic) = de(j,ic)
93                  d0(j,id) = de(j,id)
94                  d0(j,ie) = de(j,ie)
95                  d0(j,ig) = de(j,ig)
96               end do
97            end if
98            x(i) = x(i) + eps
99            call epitors2a (ipitors,de)
100            x(i) = old
101            do j = 1, 3
102               hessx(j,ia) = hessx(j,ia) + term*(de(j,ia)-d0(j,ia))
103               hessx(j,ib) = hessx(j,ib) + term*(de(j,ib)-d0(j,ib))
104               hessx(j,ic) = hessx(j,ic) + term*(de(j,ic)-d0(j,ic))
105               hessx(j,id) = hessx(j,id) + term*(de(j,id)-d0(j,id))
106               hessx(j,ie) = hessx(j,ie) + term*(de(j,ie)-d0(j,ie))
107               hessx(j,ig) = hessx(j,ig) + term*(de(j,ig)-d0(j,ig))
108            end do
109c
110c     find numerical y-components via perturbed structures
111c
112            old = y(i)
113            if (twosided) then
114               y(i) = y(i) - 0.5d0*eps
115               call epitors2a (ipitors,de)
116               do j = 1, 3
117                  d0(j,ia) = de(j,ia)
118                  d0(j,ib) = de(j,ib)
119                  d0(j,ic) = de(j,ic)
120                  d0(j,id) = de(j,id)
121                  d0(j,ie) = de(j,ie)
122                  d0(j,ig) = de(j,ig)
123               end do
124            end if
125            y(i) = y(i) + eps
126            call epitors2a (ipitors,de)
127            y(i) = old
128            do j = 1, 3
129               hessy(j,ia) = hessy(j,ia) + term*(de(j,ia)-d0(j,ia))
130               hessy(j,ib) = hessy(j,ib) + term*(de(j,ib)-d0(j,ib))
131               hessy(j,ic) = hessy(j,ic) + term*(de(j,ic)-d0(j,ic))
132               hessy(j,id) = hessy(j,id) + term*(de(j,id)-d0(j,id))
133               hessy(j,ie) = hessy(j,ie) + term*(de(j,ie)-d0(j,ie))
134               hessy(j,ig) = hessy(j,ig) + term*(de(j,ig)-d0(j,ig))
135            end do
136c
137c     find numerical z-components via perturbed structures
138c
139            old = z(i)
140            if (twosided) then
141               z(i) = z(i) - 0.5d0*eps
142               call epitors2a (ipitors,de)
143               do j = 1, 3
144                  d0(j,ia) = de(j,ia)
145                  d0(j,ib) = de(j,ib)
146                  d0(j,ic) = de(j,ic)
147                  d0(j,id) = de(j,id)
148                  d0(j,ie) = de(j,ie)
149                  d0(j,ig) = de(j,ig)
150               end do
151            end if
152            z(i) = z(i) + eps
153            call epitors2a (ipitors,de)
154            z(i) = old
155            do j = 1, 3
156               hessz(j,ia) = hessz(j,ia) + term*(de(j,ia)-d0(j,ia))
157               hessz(j,ib) = hessz(j,ib) + term*(de(j,ib)-d0(j,ib))
158               hessz(j,ic) = hessz(j,ic) + term*(de(j,ic)-d0(j,ic))
159               hessz(j,id) = hessz(j,id) + term*(de(j,id)-d0(j,id))
160               hessz(j,ie) = hessz(j,ie) + term*(de(j,ie)-d0(j,ie))
161               hessz(j,ig) = hessz(j,ig) + term*(de(j,ig)-d0(j,ig))
162            end do
163         end if
164      end do
165c
166c     perform deallocation of some local arrays
167c
168      deallocate (d0)
169      return
170      end
171c
172c
173c     ###############################################################
174c     ##                                                           ##
175c     ##  subroutine epitors2a  --  pi-system torsion derivatives  ##
176c     ##                                                           ##
177c     ###############################################################
178c
179c
180c     "epitors2a" calculates the pi-system torsion first derivatives;
181c     used in computation of finite difference second derivatives
182c
183c
184      subroutine epitors2a (i,de)
185      use atoms
186      use bound
187      use deriv
188      use pitors
189      use torpot
190      implicit none
191      integer i,ia,ib,ic
192      integer id,ie,ig
193      real*8 dedphi,dphi2
194      real*8 xt,yt,zt,rt2
195      real*8 xu,yu,zu,ru2
196      real*8 xtu,ytu,ztu
197      real*8 rdc,rtru
198      real*8 v2,c2,s2
199      real*8 sine,cosine
200      real*8 sine2,cosine2
201      real*8 xia,yia,zia
202      real*8 xib,yib,zib
203      real*8 xic,yic,zic
204      real*8 xid,yid,zid
205      real*8 xie,yie,zie
206      real*8 xig,yig,zig
207      real*8 xip,yip,zip
208      real*8 xiq,yiq,ziq
209      real*8 xad,yad,zad
210      real*8 xbd,ybd,zbd
211      real*8 xec,yec,zec
212      real*8 xgc,ygc,zgc
213      real*8 xcp,ycp,zcp
214      real*8 xdc,ydc,zdc
215      real*8 xqd,yqd,zqd
216      real*8 xdp,ydp,zdp
217      real*8 xqc,yqc,zqc
218      real*8 dedxt,dedyt,dedzt
219      real*8 dedxu,dedyu,dedzu
220      real*8 dedxia,dedyia,dedzia
221      real*8 dedxib,dedyib,dedzib
222      real*8 dedxic,dedyic,dedzic
223      real*8 dedxid,dedyid,dedzid
224      real*8 dedxie,dedyie,dedzie
225      real*8 dedxig,dedyig,dedzig
226      real*8 dedxip,dedyip,dedzip
227      real*8 dedxiq,dedyiq,dedziq
228      real*8 de(3,*)
229c
230c
231c     set the atom numbers for this pi-system torsion
232c
233      ia = ipit(1,i)
234      ib = ipit(2,i)
235      ic = ipit(3,i)
236      id = ipit(4,i)
237      ie = ipit(5,i)
238      ig = ipit(6,i)
239c
240c     compute the value of the pi-system torsion angle
241c
242      xia = x(ia)
243      yia = y(ia)
244      zia = z(ia)
245      xib = x(ib)
246      yib = y(ib)
247      zib = z(ib)
248      xic = x(ic)
249      yic = y(ic)
250      zic = z(ic)
251      xid = x(id)
252      yid = y(id)
253      zid = z(id)
254      xie = x(ie)
255      yie = y(ie)
256      zie = z(ie)
257      xig = x(ig)
258      yig = y(ig)
259      zig = z(ig)
260      xad = xia - xid
261      yad = yia - yid
262      zad = zia - zid
263      xbd = xib - xid
264      ybd = yib - yid
265      zbd = zib - zid
266      xec = xie - xic
267      yec = yie - yic
268      zec = zie - zic
269      xgc = xig - xic
270      ygc = yig - yic
271      zgc = zig - zic
272      if (use_polymer) then
273         call image (xad,yad,zad)
274         call image (xbd,ybd,zbd)
275         call image (xec,yec,zec)
276         call image (xgc,ygc,zgc)
277      end if
278      xip = yad*zbd - ybd*zad + xic
279      yip = zad*xbd - zbd*xad + yic
280      zip = xad*ybd - xbd*yad + zic
281      xiq = yec*zgc - ygc*zec + xid
282      yiq = zec*xgc - zgc*xec + yid
283      ziq = xec*ygc - xgc*yec + zid
284      xcp = xic - xip
285      ycp = yic - yip
286      zcp = zic - zip
287      xdc = xid - xic
288      ydc = yid - yic
289      zdc = zid - zic
290      xqd = xiq - xid
291      yqd = yiq - yid
292      zqd = ziq - zid
293      if (use_polymer) then
294         call image (xcp,ycp,zcp)
295         call image (xdc,ydc,zdc)
296         call image (xqd,yqd,zqd)
297      end if
298      xt = ycp*zdc - ydc*zcp
299      yt = zcp*xdc - zdc*xcp
300      zt = xcp*ydc - xdc*ycp
301      xu = ydc*zqd - yqd*zdc
302      yu = zdc*xqd - zqd*xdc
303      zu = xdc*yqd - xqd*ydc
304      xtu = yt*zu - yu*zt
305      ytu = zt*xu - zu*xt
306      ztu = xt*yu - xu*yt
307      rt2 = xt*xt + yt*yt + zt*zt
308      ru2 = xu*xu + yu*yu + zu*zu
309      rtru = sqrt(rt2 * ru2)
310      if (rtru .ne. 0.0d0) then
311         rdc = sqrt(xdc*xdc + ydc*ydc + zdc*zdc)
312         cosine = (xt*xu + yt*yu + zt*zu) / rtru
313         sine = (xdc*xtu + ydc*ytu + zdc*ztu) / (rdc*rtru)
314c
315c     set the pi-system torsion parameters for this angle
316c
317         v2 = kpit(i)
318         c2 = -1.0d0
319         s2 = 0.0d0
320c
321c     compute the multiple angle trigonometry and the phase terms
322c
323         cosine2 = cosine*cosine - sine*sine
324         sine2 = 2.0d0 * cosine * sine
325         dphi2 = 2.0d0 * (cosine2*s2 - sine2*c2)
326c
327c     calculate pi-system torsion energy and master chain rule term
328c
329         dedphi = ptorunit * v2 * dphi2
330c
331c     chain rule terms for first derivative components
332c
333         xdp = xid - xip
334         ydp = yid - yip
335         zdp = zid - zip
336         xqc = xiq - xic
337         yqc = yiq - yic
338         zqc = ziq - zic
339         dedxt = dedphi * (yt*zdc - ydc*zt) / (rt2*rdc)
340         dedyt = dedphi * (zt*xdc - zdc*xt) / (rt2*rdc)
341         dedzt = dedphi * (xt*ydc - xdc*yt) / (rt2*rdc)
342         dedxu = -dedphi * (yu*zdc - ydc*zu) / (ru2*rdc)
343         dedyu = -dedphi * (zu*xdc - zdc*xu) / (ru2*rdc)
344         dedzu = -dedphi * (xu*ydc - xdc*yu) / (ru2*rdc)
345c
346c     compute first derivative components for pi-system angle
347c
348         dedxip = zdc*dedyt - ydc*dedzt
349         dedyip = xdc*dedzt - zdc*dedxt
350         dedzip = ydc*dedxt - xdc*dedyt
351         dedxic = ydp*dedzt - zdp*dedyt + zqd*dedyu - yqd*dedzu
352         dedyic = zdp*dedxt - xdp*dedzt + xqd*dedzu - zqd*dedxu
353         dedzic = xdp*dedyt - ydp*dedxt + yqd*dedxu - xqd*dedyu
354         dedxid = zcp*dedyt - ycp*dedzt + yqc*dedzu - zqc*dedyu
355         dedyid = xcp*dedzt - zcp*dedxt + zqc*dedxu - xqc*dedzu
356         dedzid = ycp*dedxt - xcp*dedyt + xqc*dedyu - yqc*dedxu
357         dedxiq = zdc*dedyu - ydc*dedzu
358         dedyiq = xdc*dedzu - zdc*dedxu
359         dedziq = ydc*dedxu - xdc*dedyu
360c
361c     compute first derivative components for individual atoms
362c
363         dedxia = ybd*dedzip - zbd*dedyip
364         dedyia = zbd*dedxip - xbd*dedzip
365         dedzia = xbd*dedyip - ybd*dedxip
366         dedxib = zad*dedyip - yad*dedzip
367         dedyib = xad*dedzip - zad*dedxip
368         dedzib = yad*dedxip - xad*dedyip
369         dedxie = ygc*dedziq - zgc*dedyiq
370         dedyie = zgc*dedxiq - xgc*dedziq
371         dedzie = xgc*dedyiq - ygc*dedxiq
372         dedxig = zec*dedyiq - yec*dedziq
373         dedyig = xec*dedziq - zec*dedxiq
374         dedzig = yec*dedxiq - xec*dedyiq
375         dedxic = dedxic + dedxip - dedxie - dedxig
376         dedyic = dedyic + dedyip - dedyie - dedyig
377         dedzic = dedzic + dedzip - dedzie - dedzig
378         dedxid = dedxid + dedxiq - dedxia - dedxib
379         dedyid = dedyid + dedyiq - dedyia - dedyib
380         dedzid = dedzid + dedziq - dedzia - dedzib
381c
382c     increment the total pi-system torsion energy and gradient
383c
384         de(1,ia) = dedxia
385         de(2,ia) = dedyia
386         de(3,ia) = dedzia
387         de(1,ib) = dedxib
388         de(2,ib) = dedyib
389         de(3,ib) = dedzib
390         de(1,ic) = dedxic
391         de(2,ic) = dedyic
392         de(3,ic) = dedzic
393         de(1,id) = dedxid
394         de(2,id) = dedyid
395         de(3,id) = dedzid
396         de(1,ie) = dedxie
397         de(2,ie) = dedyie
398         de(3,ie) = dedzie
399         de(1,ig) = dedxig
400         de(2,ig) = dedyig
401         de(3,ig) = dedzig
402      end if
403      return
404      end
405