1c
2c
3c     ############################################################
4c     ##  COPYRIGHT (C) 1995 by Yong Kong & Jay William Ponder  ##
5c     ##                  All Rights Reserved                   ##
6c     ############################################################
7c
8c     #################################################################
9c     ##                                                             ##
10c     ##  subroutine rotpole  --  rotate multipoles to global frame  ##
11c     ##                                                             ##
12c     #################################################################
13c
14c
15c     "rotpole" constructs the set of atomic multipoles in the global
16c     frame by applying the correct rotation matrix for each site
17c
18c
19      subroutine rotpole
20      use mpole
21      implicit none
22      integer i
23      real*8 a(3,3)
24c     real*8 d1(3,3)
25c     real*8 d2(5,5)
26c     logical use_qi
27c
28c
29c     rotate the atomic multipoles at each site in turn
30c
31      do i = 1, npole
32         call rotmat (i,a)
33         call rotsite (i,a)
34      end do
35c
36c     set QI dipole rotation matrix by permuting the Cartesian
37c     rotation matrix to adhere to spherical harmonic ordering
38c
39c     do i = 1, npole
40c        d1(1,1) = a(3,3)
41c        d1(2,1) = a(3,1)
42c        d1(3,1) = a(3,2)
43c        d1(1,2) = a(1,3)
44c        d1(2,2) = a(1,1)
45c        d1(3,2) = a(1,2)
46c        d1(1,3) = a(2,3)
47c        d1(2,3) = a(2,1)
48c        d1(3,3) = a(2,2)
49c        call shrotmat (d1,d2)
50c        call shrotsite (i,d1,d2)
51c     end do
52      return
53      end
54c
55c
56c     ################################################################
57c     ##                                                            ##
58c     ##  subroutine rotmat  --  find global frame rotation matrix  ##
59c     ##                                                            ##
60c     ################################################################
61c
62c
63c     "rotmat" finds the rotation matrix that rotates the local
64c     coordinate system into the global frame at a multipole site
65c
66c
67      subroutine rotmat (i,a)
68      use atoms
69      use mpole
70      implicit none
71      integer i,ii
72      integer ix,iy,iz
73      real*8 r,dot
74      real*8 xi,yi,zi
75      real*8 dx,dy,dz
76      real*8 dx1,dy1,dz1
77      real*8 dx2,dy2,dz2
78      real*8 dx3,dy3,dz3
79      real*8 a(3,3)
80c
81c
82c     get coordinates and frame definition for the multipole site
83c
84      ii = ipole(i)
85      xi = x(ii)
86      yi = y(ii)
87      zi = z(ii)
88      iz = zaxis(i)
89      ix = xaxis(i)
90      iy = abs(yaxis(i))
91c
92c     use the identity matrix as the default rotation matrix
93c
94      a(1,1) = 1.0d0
95      a(2,1) = 0.0d0
96      a(3,1) = 0.0d0
97      a(1,3) = 0.0d0
98      a(2,3) = 0.0d0
99      a(3,3) = 1.0d0
100c
101c     z-only frame rotation matrix elements for z-axis only
102c
103      if (polaxe(i) .eq. 'Z-Only') then
104         dx = x(iz) - xi
105         dy = y(iz) - yi
106         dz = z(iz) - zi
107         r = sqrt(dx*dx + dy*dy + dz*dz)
108         a(1,3) = dx / r
109         a(2,3) = dy / r
110         a(3,3) = dz / r
111         dx = 1.0d0
112         dy = 0.0d0
113         dz = 0.0d0
114         dot = a(1,3)
115         if (abs(dot) .gt. 0.866d0) then
116            dx = 0.0d0
117            dy = 1.0d0
118            dot = a(2,3)
119         end if
120         dx = dx - dot*a(1,3)
121         dy = dy - dot*a(2,3)
122         dz = dz - dot*a(3,3)
123         r = sqrt(dx*dx + dy*dy + dz*dz)
124         a(1,1) = dx / r
125         a(2,1) = dy / r
126         a(3,1) = dz / r
127c
128c     z-then-x frame rotation matrix elements for z- and x-axes
129c
130      else if (polaxe(i) .eq. 'Z-then-X') then
131         dx = x(iz) - xi
132         dy = y(iz) - yi
133         dz = z(iz) - zi
134         r = sqrt(dx*dx + dy*dy + dz*dz)
135         a(1,3) = dx / r
136         a(2,3) = dy / r
137         a(3,3) = dz / r
138         dx = x(ix) - xi
139         dy = y(ix) - yi
140         dz = z(ix) - zi
141         dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3)
142         dx = dx - dot*a(1,3)
143         dy = dy - dot*a(2,3)
144         dz = dz - dot*a(3,3)
145         r = sqrt(dx*dx + dy*dy + dz*dz)
146         a(1,1) = dx / r
147         a(2,1) = dy / r
148         a(3,1) = dz / r
149c
150c     bisector frame rotation matrix elements for z- and x-axes
151c
152      else if (polaxe(i) .eq. 'Bisector') then
153         dx = x(iz) - xi
154         dy = y(iz) - yi
155         dz = z(iz) - zi
156         r = sqrt(dx*dx + dy*dy + dz*dz)
157         dx1 = dx / r
158         dy1 = dy / r
159         dz1 = dz / r
160         dx = x(ix) - xi
161         dy = y(ix) - yi
162         dz = z(ix) - zi
163         r = sqrt(dx*dx + dy*dy + dz*dz)
164         dx2 = dx / r
165         dy2 = dy / r
166         dz2 = dz / r
167         dx = dx1 + dx2
168         dy = dy1 + dy2
169         dz = dz1 + dz2
170         r = sqrt(dx*dx + dy*dy + dz*dz)
171         a(1,3) = dx / r
172         a(2,3) = dy / r
173         a(3,3) = dz / r
174         dot = dx2*a(1,3) + dy2*a(2,3) + dz2*a(3,3)
175         dx = dx2 - dot*a(1,3)
176         dy = dy2 - dot*a(2,3)
177         dz = dz2 - dot*a(3,3)
178         r = sqrt(dx*dx + dy*dy + dz*dz)
179         a(1,1) = dx / r
180         a(2,1) = dy / r
181         a(3,1) = dz / r
182c
183c     z-bisect frame rotation matrix elements for z- and x-axes
184c
185      else if (polaxe(i) .eq. 'Z-Bisect') then
186         dx = x(iz) - xi
187         dy = y(iz) - yi
188         dz = z(iz) - zi
189         r = sqrt(dx*dx + dy*dy + dz*dz)
190         a(1,3) = dx / r
191         a(2,3) = dy / r
192         a(3,3) = dz / r
193         dx = x(ix) - xi
194         dy = y(ix) - yi
195         dz = z(ix) - zi
196         r = sqrt(dx*dx + dy*dy + dz*dz)
197         dx1 = dx / r
198         dy1 = dy / r
199         dz1 = dz / r
200         dx = x(iy) - xi
201         dy = y(iy) - yi
202         dz = z(iy) - zi
203         r = sqrt(dx*dx + dy*dy + dz*dz)
204         dx2 = dx / r
205         dy2 = dy / r
206         dz2 = dz / r
207         dx = dx1 + dx2
208         dy = dy1 + dy2
209         dz = dz1 + dz2
210         r = sqrt(dx*dx + dy*dy + dz*dz)
211         dx = dx / r
212         dy = dy / r
213         dz = dz / r
214         dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3)
215         dx = dx - dot*a(1,3)
216         dy = dy - dot*a(2,3)
217         dz = dz - dot*a(3,3)
218         r = sqrt(dx*dx + dy*dy + dz*dz)
219         a(1,1) = dx / r
220         a(2,1) = dy / r
221         a(3,1) = dz / r
222c
223c     3-fold frame rotation matrix elements for z- and x-axes
224c
225      else if (polaxe(i) .eq. '3-Fold') then
226         dx = x(iz) - xi
227         dy = y(iz) - yi
228         dz = z(iz) - zi
229         r = sqrt(dx*dx + dy*dy + dz*dz)
230         dx1 = dx / r
231         dy1 = dy / r
232         dz1 = dz / r
233         dx = x(ix) - xi
234         dy = y(ix) - yi
235         dz = z(ix) - zi
236         r = sqrt(dx*dx + dy*dy + dz*dz)
237         dx2 = dx / r
238         dy2 = dy / r
239         dz2 = dz / r
240         dx = x(iy) - xi
241         dy = y(iy) - yi
242         dz = z(iy) - zi
243         r = sqrt(dx*dx + dy*dy + dz*dz)
244         dx3 = dx / r
245         dy3 = dy / r
246         dz3 = dz / r
247         dx = dx1 + dx2 + dx3
248         dy = dy1 + dy2 + dy3
249         dz = dz1 + dz2 + dz3
250         r = sqrt(dx*dx + dy*dy + dz*dz)
251         a(1,3) = dx / r
252         a(2,3) = dy / r
253         a(3,3) = dz / r
254         dot = dx2*a(1,3) + dy2*a(2,3) + dz2*a(3,3)
255         dx = dx2 - dot*a(1,3)
256         dy = dy2 - dot*a(2,3)
257         dz = dz2 - dot*a(3,3)
258         r = sqrt(dx*dx + dy*dy + dz*dz)
259         a(1,1) = dx / r
260         a(2,1) = dy / r
261         a(3,1) = dz / r
262      end if
263c
264c     finally, find rotation matrix elements for the y-axis
265c
266      a(1,2) = a(3,1)*a(2,3) - a(2,1)*a(3,3)
267      a(2,2) = a(1,1)*a(3,3) - a(3,1)*a(1,3)
268      a(3,2) = a(2,1)*a(1,3) - a(1,1)*a(2,3)
269      return
270      end
271c
272c
273c     ################################################################
274c     ##                                                            ##
275c     ##  subroutine rotsite  --  rotate multipoles at single site  ##
276c     ##                                                            ##
277c     ################################################################
278c
279c
280c     "rotsite" rotates the local frame atomic multipoles at a
281c     specified site into the global coordinate frame by applying
282c     a rotation matrix
283c
284c
285      subroutine rotsite (isite,a)
286      use atoms
287      use mpole
288      implicit none
289      integer i,j,k,m
290      integer isite
291      real*8 a(3,3)
292      real*8 mp(3,3)
293      real*8 rp(3,3)
294c
295c
296c     monopoles have the same value in any coordinate frame
297c
298      rpole(1,isite) = pole(1,isite)
299c
300c     rotate the dipoles to the global coordinate frame
301c
302      do i = 2, 4
303         rpole(i,isite) = 0.0d0
304         do j = 2, 4
305            rpole(i,isite) = rpole(i,isite) + pole(j,isite)*a(i-1,j-1)
306         end do
307      end do
308c
309c     rotate the quadrupoles to the global coordinate frame
310c
311      k = 5
312      do i = 1, 3
313         do j = 1, 3
314            mp(i,j) = pole(k,isite)
315            rp(i,j) = 0.0d0
316            k = k + 1
317         end do
318      end do
319      do i = 1, 3
320         do j = 1, 3
321            if (j .lt. i) then
322               rp(i,j) = rp(j,i)
323            else
324               do k = 1, 3
325                  do m = 1, 3
326                     rp(i,j) = rp(i,j) + a(i,k)*a(j,m)*mp(k,m)
327                  end do
328               end do
329            end if
330         end do
331      end do
332      k = 5
333      do i = 1, 3
334         do j = 1, 3
335            rpole(k,isite) = rp(i,j)
336            k = k + 1
337         end do
338      end do
339      return
340      end
341c
342c
343c     ##############################################################
344c     ##                                                          ##
345c     ##  subroutine shrotmat  --  QI quadrupole rotation matrix  ##
346c     ##                                                          ##
347c     ##############################################################
348c
349c
350c     "shrotmat" finds the rotation matrix that converts spherical
351c     harmonic quadrupoles from the local to the global frame given
352c     the required dipole rotation matrix
353c
354c
355      subroutine shrotmat (d1,d2)
356      use math
357      implicit none
358      real*8 d1(3,3)
359      real*8 d2(5,5)
360c
361c
362c     build the quadrupole rotation matrix
363c
364      d2(1,1) = 0.5d0*(3.0d0*d1(1,1)*d1(1,1) - 1.0d0)
365      d2(2,1) = root3*d1(1,1)*d1(2,1)
366      d2(3,1) = root3*d1(1,1)*d1(3,1)
367      d2(4,1) = 0.5d0*root3*(d1(2,1)*d1(2,1) - d1(3,1)*d1(3,1))
368      d2(5,1) = root3*d1(2,1)*d1(3,1)
369      d2(1,2) = root3*d1(1,1)*d1(1,2)
370      d2(2,2) = d1(2,1)*d1(1,2) + d1(1,1)*d1(2,2)
371      d2(3,2) = d1(3,1)*d1(1,2) + d1(1,1)*d1(3,2)
372      d2(4,2) = d1(2,1)*d1(2,2) - d1(3,1)*d1(3,2)
373      d2(5,2) = d1(3,1)*d1(2,2) + d1(2,1)*d1(3,2)
374      d2(1,3) = root3*d1(1,1)*d1(1,3)
375      d2(2,3) = d1(2,1)*d1(1,3) + d1(1,1)*d1(2,3)
376      d2(3,3) = d1(3,1)*d1(1,3) + d1(1,1)*d1(3,3)
377      d2(4,3) = d1(2,1)*d1(2,3) - d1(3,1)*d1(3,3)
378      d2(5,3) = d1(3,1)*d1(2,3) + d1(2,1)*d1(3,3)
379      d2(1,4) = 0.5d0*root3*(d1(1,2)*d1(1,2) - d1(1,3)*d1(1,3))
380      d2(2,4) = d1(1,2)*d1(2,2) - d1(1,3)*d1(2,3)
381      d2(3,4) = d1(1,2)*d1(3,2) - d1(1,3)*d1(3,3)
382      d2(4,4) = 0.5d0*(d1(2,2)*d1(2,2) - d1(3,2)*d1(3,2)
383     &             - d1(2,3)*d1(2,3) + d1(3,3)*d1(3,3))
384      d2(5,4) = d1(2,2)*d1(3,2) - d1(2,3)*d1(3,3)
385      d2(1,5) = root3*d1(1,2)*d1(1,3)
386      d2(2,5) = d1(2,2)*d1(1,3) + d1(1,2)*d1(2,3)
387      d2(3,5) = d1(3,2)*d1(1,3) + d1(1,2)*d1(3,3)
388      d2(4,5) = d1(2,2)*d1(2,3) - d1(3,2)*d1(3,3)
389      d2(5,5) = d1(3,2)*d1(2,3) + d1(2,2)*d1(3,3)
390      return
391      end
392c
393c
394c     ##################################################################
395c     ##                                                              ##
396c     ##  subroutine shrotsite  --  rotate SH mpoles local to global  ##
397c     ##                                                              ##
398c     ##################################################################
399c
400c
401c     "shrotsite" converts spherical harmonic multipoles from the
402c     local to the global frame given required rotation matrices
403c
404c
405      subroutine shrotsite (i,d1,d2)
406      use mpole
407      implicit none
408      integer i
409      real*8 d1(3,3)
410      real*8 d2(5,5)
411c
412c
413c     find the global frame spherical harmonic multipoles
414c
415      srpole(1,i) = spole(1,i)
416      srpole(2,i) = spole(2,i)*d1(1,1) + spole(3,i)*d1(2,1)
417     &                 + spole(4,i)*d1(3,1)
418      srpole(3,i) = spole(2,i)*d1(1,2) + spole(3,i)*d1(2,2)
419     &                 + spole(4,i)*d1(3,2)
420      srpole(4,i) = spole(2,i)*d1(1,3) + spole(3,i)*d1(2,3)
421     &                 + spole(4,i)*d1(3,3)
422      srpole(5,i) = spole(5,i)*d2(1,1) + spole(6,i)*d2(2,1)
423     &                 + spole(7,i)*d2(3,1) + spole(8,i)*d2(4,1)
424     &                 + spole(9,i)*d2(5,1)
425      srpole(6,i) = spole(5,i)*d2(1,2) + spole(6,i)*d2(2,2)
426     &                 + spole(7,i)*d2(3,2) + spole(8,i)*d2(4,2)
427     &                 + spole(9,i)*d2(5,2)
428      srpole(7,i) = spole(5,i)*d2(1,3) + spole(6,i)*d2(2,3)
429     &                 + spole(7,i)*d2(3,3) + spole(8,i)*d2(4,3)
430     &                 + spole(9,i)*d2(5,3)
431      srpole(8,i) = spole(5,i)*d2(1,4) + spole(6,i)*d2(2,4)
432     &                 + spole(7,i)*d2(3,4) + spole(8,i)*d2(4,4)
433     &                 + spole(9,i)*d2(5,4)
434      srpole(9,i) = spole(5,i)*d2(1,5) + spole(6,i)*d2(2,5)
435     &                 + spole(7,i)*d2(3,5) + spole(8,i)*d2(4,5)
436     &                 + spole(9,i)*d2(5,5)
437      return
438      end
439c
440c
441c     ###############################################################
442c     ##                                                           ##
443c     ##  subroutine qirotmat  --  QI interatomic rotation matrix  ##
444c     ##                                                           ##
445c     ###############################################################
446c
447c
448c     "qirotmat" finds a rotation matrix that describes the
449c     interatomic vector
450c
451c
452      subroutine qirotmat (i,k,rinv,d1)
453      use atoms
454      use mpole
455      implicit none
456      integer i,k
457      real*8 rinv,r,dot
458      real*8 dx,dy,dz
459      real*8 a(3,3)
460      real*8 d1(3,3)
461c
462c     set the z-axis to be the interatomic vector
463c
464      dx = x(k) - x(i)
465      dy = y(k) - y(i)
466      dz = z(k) - z(i)
467      a(1,3) = dx * rinv
468      a(2,3) = dy * rinv
469      a(3,3) = dz * rinv
470c
471c     find an x-axis that is orthogonal to the z-axis
472c
473      if (y(i).ne.y(k) .or. z(i).ne.z(k)) then
474        dx = dx + 1.0d0
475      else
476        dy = dy + 1.0d0
477      end if
478      dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3)
479      dx = dx - dot*a(1,3)
480      dy = dy - dot*a(2,3)
481      dz = dz - dot*a(3,3)
482      r = 1.d0 / sqrt(dx*dx + dy*dy + dz*dz)
483      a(1,1) = dx * r
484      a(2,1) = dy * r
485      a(3,1) = dz * r
486c
487c     get rotation matrix elements for the y-axis
488c
489      a(1,2) = a(3,1)*a(2,3) - a(2,1)*a(3,3)
490      a(2,2) = a(1,1)*a(3,3) - a(3,1)*a(1,3)
491      a(3,2) = a(2,1)*a(1,3) - a(1,1)*a(2,3)
492c
493c     reorder to account for spherical harmonics
494c
495      d1(1,1) = a(3,3)
496      d1(2,1) = a(1,3)
497      d1(3,1) = a(2,3)
498      d1(1,2) = a(3,1)
499      d1(2,2) = a(1,1)
500      d1(3,2) = a(2,1)
501      d1(1,3) = a(3,2)
502      d1(2,3) = a(1,2)
503      d1(3,3) = a(2,2)
504      return
505      end
506