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