1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1992  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ##################################################################
9c     ##                                                              ##
10c     ##  subroutine switch  --  get switching function coefficients  ##
11c     ##                                                              ##
12c     ##################################################################
13c
14c
15c     "switch" sets the coeffcients used by the fifth and seventh
16c     order polynomial switching functions for spherical cutoffs
17c
18c
19      subroutine switch (mode)
20      use limits
21      use nonpol
22      use shunt
23      implicit none
24      real*8 denom,term
25      real*8 off3,off4,off5
26      real*8 off6,off7
27      real*8 cut3,cut4,cut5
28      real*8 cut6,cut7
29      character*6 mode
30c
31c
32c     get the switching window for the current potential type
33c
34      if (mode(1:3) .eq. 'VDW') then
35         off = vdwcut
36         cut = vdwtaper
37      else if (mode(1:6) .eq. 'REPULS') then
38         off = repcut
39         cut = reptaper
40      else if (mode(1:4) .eq. 'DISP') then
41         off = dispcut
42         cut = disptaper
43      else if (mode(1:6) .eq. 'CHARGE') then
44         off = chgcut
45         cut = chgtaper
46      else if (mode(1:6) .eq. 'CHGDPL') then
47         off = sqrt(chgcut*dplcut)
48         cut = sqrt(chgtaper*dpltaper)
49      else if (mode(1:6) .eq. 'DIPOLE') then
50         off = dplcut
51         cut = dpltaper
52      else if (mode(1:5) .eq. 'MPOLE') then
53         off = mpolecut
54         cut = mpoletaper
55      else if (mode(1:6) .eq. 'CHGTRN') then
56         off = ctrncut
57         cut = ctrntaper
58      else if (mode(1:5) .eq. 'EWALD') then
59         off = ewaldcut
60         cut = ewaldcut
61      else if (mode(1:6) .eq. 'DEWALD') then
62         off = dewaldcut
63         cut = dewaldcut
64      else if (mode(1:6) .eq. 'USOLVE') then
65         off = usolvcut
66         cut = usolvcut
67      else if (mode(1:3) .eq. 'GKV') then
68         off = spoff
69         cut = spcut
70      else if (mode(1:4) .eq. 'GKSA') then
71         off = stcut
72         cut = stoff
73      else
74         off = min(vdwcut,repcut,dispcut,chgcut,
75     &                dplcut,mpolecut,ctrncut)
76         cut = min(vdwtaper,reptaper,disptaper,chgtaper,
77     &                dpltaper,mpoletaper,ctrntaper)
78      end if
79c
80c     test for replicate periodic boundaries at this cutoff
81c
82      call replica (off)
83c
84c     set switching coefficients to zero for truncation cutoffs
85c
86      c0 = 0.0d0
87      c1 = 0.0d0
88      c2 = 0.0d0
89      c3 = 0.0d0
90      c4 = 0.0d0
91      c5 = 0.0d0
92      f0 = 0.0d0
93      f1 = 0.0d0
94      f2 = 0.0d0
95      f3 = 0.0d0
96      f4 = 0.0d0
97      f5 = 0.0d0
98      f6 = 0.0d0
99      f7 = 0.0d0
100c
101c     store the powers of the switching window cutoffs
102c
103      off2 = off * off
104      off3 = off2 * off
105      off4 = off2 * off2
106      off5 = off2 * off3
107      off6 = off3 * off3
108      off7 = off3 * off4
109      cut2 = cut * cut
110      cut3 = cut2 * cut
111      cut4 = cut2 * cut2
112      cut5 = cut2 * cut3
113      cut6 = cut3 * cut3
114      cut7 = cut3 * cut4
115c
116c     get 5th degree multiplicative switching function coefficients
117c
118      if (cut .lt. off) then
119         denom = (off-cut)**5
120         c0 = off*off2 * (off2-5.0d0*off*cut+10.0d0*cut2) / denom
121         c1 = -30.0d0 * off2*cut2 / denom
122         c2 = 30.0d0 * (off2*cut+off*cut2) / denom
123         c3 = -10.0d0 * (off2+4.0d0*off*cut+cut2) / denom
124         c4 = 15.0d0 * (off+cut) / denom
125         c5 = -6.0d0 / denom
126      end if
127c
128c     get 7th degree additive switching function coefficients
129c
130      if (cut.lt.off .and. mode(1:6).eq.'CHARGE') then
131         term = 9.3d0 * cut*off / (off-cut)
132         denom = cut7 - 7.0d0*cut6*off + 21.0d0*cut5*off2
133     &              - 35.0d0*cut4*off3 + 35.0d0*cut3*off4
134     &              - 21.0d0*cut2*off5 + 7.0d0*cut*off6 - off7
135         denom = term * denom
136         f0 = cut3*off3 * (-39.0d0*cut+64.0d0*off) / denom
137         f1 = cut2*off2
138     &           * (117.0d0*cut2-100.0d0*cut*off-192.0d0*off2) / denom
139         f2 = cut*off * (-117.0d0*cut3-84.0d0*cut2*off
140     &                   +534.0d0*cut*off2+192.0d0*off3) / denom
141         f3 = (39.0d0*cut4+212.0d0*cut3*off-450.0d0*cut2*off2
142     &            -612.0d0*cut*off3-64.0d0*off4) / denom
143         f4 = (-92.0d0*cut3+66.0d0*cut2*off
144     &            +684.0d0*cut*off2+217.0d0*off3) / denom
145         f5 = (42.0d0*cut2-300.0d0*cut*off-267.0d0*off2) / denom
146         f6 = (36.0d0*cut+139.0d0*off) / denom
147         f7 = -25.0d0 / denom
148      end if
149      return
150      end
151