C==================================================== SUBROUTINE EDANI (EDA, WHICH) C C Calls EDANI2, which does the actual energy calculation C C by Nico Tjandra June 1997 C==================================================== IMPLICIT NONE C include files INCLUDE 'cns.inc' INCLUDE 'diff_anis.inc' INCLUDE 'heap.inc' C i/o DOUBLE PRECISION EDA CHARACTER*7 WHICH C IF (DANIIPTR.NE.0) THEN CALL EDANI2(EDA, HEAP(DANIIPTR), HEAP(DANIJPTR), HEAP(DANIKPTR), & HEAP(DANILPTR), HEAP(DANIMPTR), HEAP(DANINPTR), & HEAP(DANIOBSPTR), HEAP(DANIERRPTR), & HEAP(CALCDANIPTR), WHICH) END IF RETURN END C C===================================================== SUBROUTINE EDANI2 (EDA, ATOMI, ATOMJ, ATOMK, ATOML, & ATOMM, ATOMN, DANIOBS, DANIERR, DANICALC, WHICH) C C Calculates anisotropy energies C C Anis energies are of the form C E = k1*deltaAnis**2 C where C k1 = force constant, C deltaAnis = calculated Anis - observed Anis C C which is a flag that switches between energy & force and C calculated Anis (for violations) calcs C C by Nico Tjandra June 1997 C C========================================================= IMPLICIT NONE C include files INCLUDE 'cns.inc' INCLUDE 'coord.inc' INCLUDE 'numbers.inc' INCLUDE 'deriv.inc' INCLUDE 'diff_anis.inc' INCLUDE 'consta.inc' INCLUDE 'mtf.inc' C i/o INTEGER ATOMI(*), ATOMJ(*), ATOMK(*), ATOML(*) INTEGER ATOMM(*), ATOMN(*) DOUBLE PRECISION DANIOBS(*), DANIERR(*), DANICALC(*) DOUBLE PRECISION EDA CHARACTER*7 WHICH C local variables INTEGER COUNT, CLASS DOUBLE PRECISION XI, XJ, XK, XL, XM, XN, YI, YJ, YK, YL, YM, YN, & ZI, ZJ, ZK, ZL, ZM, ZN, & CALCDANI, DELTADANI, E, DF1, DF2, & DT2XI, DT2YI, DT2ZI, DT2XJ, DT2YJ, DT2ZJ, & DT2XK, DT2YK, DT2ZK, DT2XL, DT2YL, DT2ZL, & DT2XM, DT2YM, DT2ZM, DT2XN, DT2YN, DT2ZN, & DT1XI, DT1YI, DT1ZI, DT1XJ, DT1YJ, DT1ZJ, & DT1XK, DT1YK, DT1ZK, DT1XL, DT1YL, DT1ZL, & DT1XM, DT1YM, DT1ZM, DT1XN, DT1YN, DT1ZN, & o1, o2, o3, o4 ,o5 ,o6 ,o7 ,o8, o9, o10, & o11, o12, o13, o14, o15, o16, o17, o18, o19, & o20, o21, o22, o23, o24, o25, o26, o27, o28, & o29, o30, o31, o32, o33, o34, o35, o36, o37, & o38, o39, o40, o41, o42, o43, o44, o45, o46, & o47, o48, o49, o50, o51, o52, o53, o54, o55, & o56, o57, o58, o59, o60, o61, o62, o63, o64, & o65, o66, o67, o68, o69, o70, o71, o72, o73 DOUBLE PRECISION o74, o75, o76, o77, o78, o79, o80, o81, o82 DOUBLE PRECISION o83, o84, o85, o86, o87, o88, o89, o90, o91, & o92, o93, o94, o95, o96, o97, o98, o99, o100, & o101, o102, o103, o104, o105, o106, o107, o108, & o109, o110, o111, o112, o113, o114, o115, o116, & o117, o118, o119, o120, o121, o122, o123, o124, & o125, o126, o127, o128, o129, o130, o131, o132, & o133, o134, o135, o136, o137, o138, o139, o140, & o141, o142, o143, o144, o145, o146, o147, o148, & o149, o150, o151, o152, o153, o154, o155, o156, & o157, o158, o159, o160, o161, o162, o163, o164, & o165, o166, o167, o168, o169, o170, o171, o172, & o173, o174, o175, o176, o177, o178, o179, o180, & o181, o182, o183, o184, o185, o186, o187, o188, & o189, o190, o191, o192, o193, o194, o195, o196, & o197, o198, o199, o200, o201, o202, o203, o204, & o205, o206, o207, o208, o209, o210, o211, o212, & o213, o214, o215, o216, o217, o218, o219, o220 DOUBLE PRECISION o221, o222, o223, o224, o225, o226, o227, o228, & o229, o230, o231, o232, o233, o234, o235, o236, & o237, o238, o239, o240, o241, o242, o243, o244, & o245, o246, o247, o248, o249, o250, o251, o252, & o253, o254, o255, o256, o257, o258, o259, o260, & o261, o262, o263, o264, o265, o266, o267, o268, & o269, o270, o271, o272, o273, o274, o275, o276, & o277, o278, o279, o280, o281, o282, o283, o284, & o285, o286, o287, o288, o289, o290, o291, o292, & o293, o294, o295, o296, o297, o298, o299, o300, & o301, o302, o303, o304, o305, o306, o307, o308, & o309, o310, o311, o312, o313, o314, o315, o316, & o317, o318, o319, o320, o321, o322, o323, o324, & o325, o326, o327, o328, o329, o330, o331, o332, & o333, o334, o335, o336, o337, o338, o339, o340, & o341, o342, o343, o344, o345, o346, o347, o348 DOUBLE PRECISION o349, o350, o351, o352, o353, o354, o355, o356, & o357, o358, o359, o360, o361, o362, o363, o364, & o365, o366, o367, o368, o369, o370, o371, o372, & o373, o374, o375, o376, o377, o378, o379, o380, & o381, o382, o383, o384, o385, o386, o387, o388, & o389, o390, o391, o392, o393, o394, o395, o396, & o397, o398, o399, o400, o401, o402, o403, o404, & o405, o406, o407, o408, o409, o410, o411, o412, & o413, o414, o415, o416, o417, o418, o419, o420, & o421, o422, o423, o424, o425, o426, o427, o428, & o429, o430, o431, o432, o433, o434, o435, o436, & o437, o438, o439, o440, o441, o442, o443, o444, & o445, o446, o447, o448, o449, o450, o451, o452, & o453, o454, o455, o456, o457, o458, o459, o460, & o461, o462, o463, o464, o465, o466, o467, o468, & o469, o470, o471, o472, o473, o474, o475, o476 DOUBLE PRECISION o477, o478, o479, o480, o481, o482, o483, o484, & o485, o486, o487, o488, o489, o490, o491, o492, & o493, o494, o495, o496, o497, o498, o499, o500, & o501, o502, o503, o504, o505, o506, o507, o508, & o509, o510, o511, o512, o513, o514, o515, o516, & o517, o518, o519, o520, o521, o522, o523, o524, & o525, o526, o527, o528, o529, o530, o531, o532, & o533, o534, o535, o536, o537, o538, o539, o540, & o541, o542, o543, o544, o545, o546, o547, o548, & o549, o550, o551, o552, o553, o554, o555, o556, & o557, o558, o559, o560, o561, o562, o563, o564, & o565, o566, o567, o568, o569, o570, o571, o572, & o573, o574, o575, o576, o577, o578, o579, o580, & o581, o582, o583, o584, o585, o586, o587, o588, & o589, o590, o591, o592, o593, o594, o595, o596, & o597, o598, o599, o600, o601, o602, o603, o604 DOUBLE PRECISION o605, o606, o607, o608, o609, o610, o611, o612, & o613, o614, o615, o616, o617, o618, o619, o620, & o621, o622, o623, o624, o625, o626, o627, o628, & o629, o630, o631, o632, o633, o634, o635, o636, & o637, o638, o639, o640, o641, o642, o643, o644, & o645, o646, o647, o648, o649, o650, o651, o652, & o653, o654, o655, o656, o657, o658, o659, o660, & o661, o662, o663, o664, o665, o666, o667, o668, & o669, o670, o671, o672, o673, o674, o675, o676, & o677, o678, o679, o680, o681, o682, o683, o684, & o685, o686, o687, o688, o689, o690, o691, o692, & o693, o694, o695, o696, o697, o698, o699, o700, & o701, o702, o703, o704, o705, o706, o707, o708, & o709, o710, o711, o712, o713, o714, o715, o716, & o717, o718, o719, o720, o721, o722, o723, o724, & o725, o726, o727, o728, o729, o730, o731, o732 DOUBLE PRECISION o733, o734, o735, o736, o737, o738, o739, o740, & o741, o742, o743, o744, o745, o746, o747, o748, & o749, o750, o751, o752, o753, o754, o755, o756, & o757, o758, o759, o760, o761, o762, o763, o764, & o765, o766, o767, o768, o769, o770, o771, o772, & o773, o774, o775, o776, o777, o778, o779, o780, & o781, o782, o783, o784, o785, o786, o787, o788, & o789, o790, o791, o792, o793, o794, o795, o796, & o797, o798, o799, o800, o801, o802, o803, o804, & o805, o806, o807, o808, o809, o810, o811, o812, & o813, o814, o815, o816, o817, o818, o819, o820, & o821, o822, o823, o824, o825, o826, o827, o828, & o829, o830, o831, o832, o833, o834, o835, o836, & o837, o838, o839, o840, o841, o842, o843, o844, & o845, o846, o847, o848, o849, o850, o851, o852, & o853, o854, o855, o856, o857, o858, o859, o860 DOUBLE PRECISION o861, o862, o863, o864, o865, o866, o867, o868, & o869, o870, o871, o872, o873, o874, o875, o876, & o877, o878, o879, o880, o881, o882, o883, o884, & o885, o886, o887, o888, o889, o890, o891, o892, & o893, o894, o895, o896, o897, o898, o899, o900, & o901, o902, o903, o904, o905, o906, o907, o908, & o909, o910, o911, o912, o913, o914, o915, o916, & o917, o918, o919, o920, o921, o922, o923, o924, & o925, o926, o927, o928, o929, o930, o931, o932, & o933, o934, o935, o936, o937, o938, o939, o940, & o941, o942, o943, o944, o945, o946, o947, o948, & o949, o950, o951, o952, o953, o954, o955, o956, & o957, o958, o959, o960, o961, o962, o963, o964, & o965, o966, o967, o968, o969, o970, o971, o972, & o973, o974, o975, o976, o977, o978, o979, o980, & o981, o982, o983, o984, o985, o986, o987, o988 DOUBLE PRECISION o989, o990, o991, o992, o993, o994, o995, o996, & o997, o998, o999, o1000, o1001, o1002, o1003, & o1004, o1005, o1006, o1007, o1008, o1009, o1010, & o1011, o1012, o1013, o1014, o1015, o1016, o1017, & o1018, o1019, o1020, o1021, o1022, o1023, o1024, & o1025, o1026, o1027, o1028, o1029, o1030, o1031, & o1032, o1033, o1034, o1035, o1036, o1037, o1038, & o1039, o1040, o1041, o1042, o1043, o1044, o1045, & o1046, o1047, o1048, o1049, o1050, o1051, o1052, & o1053, o1054, o1055, o1056, o1057, o1058, o1059, & o1060, o1061, o1062, o1063, o1064, o1065, o1066, & o1067, o1068, o1069, o1070, o1071, o1072, o1073, & o1074, o1075, o1076, o1077, o1078, o1079, o1080, & o1081, o1082, o1083, o1084, o1085, o1086, o1087, & o1088, o1089, o1090, o1091, o1092, o1093, o1094, & o1095, o1096, o1097, o1098, o1099, o1100 DOUBLE PRECISION o1101, o1102, o1103, o1104, o1105, o1106, o1107, & o1108, o1109, o1110, o1111, o1112, o1113, o1114, & o1115, o1116, o1117, o1118, o1119, o1120, o1121, & o1122, o1123, o1124, o1125, o1126, o1127, o1128, & o1129, o1130, o1131, o1132, o1133, o1134, o1135, & o1136, o1137, o1138, o1139, o1140, o1141, o1142, & o1143, o1144, o1145, o1146, o1147, o1148, o1149, & o1150, o1151, o1152, o1153, o1154, o1155, o1156, & o1157, o1158, o1159, o1160, o1161, o1162, o1163, & o1164, o1165, o1166, o1167, o1168, o1169, o1170, & o1171, o1172, o1173, o1174, o1175, o1176, o1177, & o1178, o1179, o1180, o1181, o1182, o1183, o1184, & o1185, o1186, o1187, o1188, o1189, o1190, o1191, & o1192, o1193, o1194, o1195, o1196, o1197, o1198, & o1199, o1200 DOUBLE PRECISION OBSDANI, ERRDANI, K1, COEF0, & COEF1, COEF2, COEF3, COEF4, As, An, TC, wh, wn, & h, Gh, Gn, dist, CSA, cc, dd, T1, T2, mu C C begin C C following Axel's code in ETOR, C C zero out partial energy C EDA = ZERO C CLASS = 1 K1 = DANIFORCES(1) COEF0 = DANICOEF0(1) COEF1 = DANICOEF1(1) COEF2 = DANICOEF2(1) COEF3 = DANICOEF3(1) COEF4 = DANICOEF4(1) C COUNT = 0 C DO WHILE(COUNT.LT.NDANI) COUNT = COUNT+1 C IF (DANIASSNDX(CLASS).LT.COUNT) THEN CLASS = CLASS + 1 IF (DANIASSNDX(CLASS).EQ.DANIASSNDX(CLASS-1)) THEN COUNT = COUNT-1 ELSE K1 = DANIFORCES(CLASS) COEF0 = DANICOEF0(CLASS) COEF1 = DANICOEF1(CLASS) COEF2 = DANICOEF2(CLASS) COEF3 = DANICOEF3(CLASS) COEF4 = DANICOEF4(CLASS) END IF END IF C XI = X(ATOMI(COUNT)) XJ = X(ATOMJ(COUNT)) XK = X(ATOMK(COUNT)) XL = X(ATOML(COUNT)) XM = X(ATOMM(COUNT)) XN = X(ATOMN(COUNT)) YI = Y(ATOMI(COUNT)) YJ = Y(ATOMJ(COUNT)) YK = Y(ATOMK(COUNT)) YL = Y(ATOML(COUNT)) YM = Y(ATOMM(COUNT)) YN = Y(ATOMN(COUNT)) ZI = Z(ATOMI(COUNT)) ZJ = Z(ATOMJ(COUNT)) ZK = Z(ATOMK(COUNT)) ZL = Z(ATOML(COUNT)) ZM = Z(ATOMM(COUNT)) ZN = Z(ATOMN(COUNT)) C OBSDANI = DANIOBS(COUNT) ERRDANI = DANIERR(COUNT) C C C C Define constants C h = 1.054589E-34 Gh = 2.6753E8 Gn = -2.712E7 dist = 1.02E-10 CSA = -160.E-6 mu = 1.E-7 C C C C Calculate T1 & T2 and its derivatives with respect to x,y,z C of atoms: i, j, k, l, m, and n. This was done with Mathematica !!! C C COEF0 = Correlation time Tc (ns) C COEF1 = Dz / 0.5*(Dx+Dy) C COEF2 = 1.5*(Dy-Dx)/(Dz-0.5*(Dy+Dx)) C COEF3 = Wh (MHz) C COEF4 = Wn (MHz) C C C TC = COEF0*1.E-9 An = COEF1 C As = (3.+0.5*(An*COEF2))/(3.-0.5*(An*COEF2)) As = (3.+COEF2*(An-1.))/(3.-COEF2*(An-1.)) wh = COEF3*1.E6*2.0*3.141592654 wn = -1.0*COEF4*1.E6*2.0*3.141592654 C cc = (2./15.)*((wn*CSA)**2) dd = 0.1*(Gh*Gn*h*mu/(dist**3))**2 C C C o1=An*(1.d0+As) o2=As+5.d-1*o1 o3=1/(1.d0+o2) o4=1/TC o5=o3*o4 o6=As*o3*o4 o7=An*(1.d0+As)*o3*o4 o8=5.d-1*o5+2.d0*o6+2.5d-1*o7 o9=1/o8 o10=o8**2 o11=1/o10 o12=wn**2 o13=o11*o12 o14=1/(1.d0+1.d0*o13) o15=-XI+XJ o16=o15**2 o17=-YI+YJ o18=o17**2 o19=-ZI+ZJ o20=o19**2 o21=o16+o18+o20 o22=1/o21 o23=-XI+XK o24=o23**2 o25=-YI+YK o26=o25**2 o27=-ZI+ZK o28=o27**2 o29=o24+o26+o28 o30=1/o29 o31=-XM+XN o32=o15*o31 o33=-YM+YN o34=o17*o33 o35=-ZM+ZN o36=o19*o35 o37=o32+o34+o36 o38=o37**2 o39=o23*o31 o40=o25*o33 o41=o27*o35 o42=o39+o40+o41 o43=o42**2 o44=o31**2 o45=o33**2 o46=o35**2 o47=o44+o45+o46 o48=o47**2 o49=1/o48 o50=o14*o22*o30*o38*o43*o49*o9 o51=2.d0*o5+5.d-1*o6+2.5d-1*o7 o52=1/o51 o53=o51**2 o54=1/o53 o55=o12*o54 o56=1/(1.d0+1.d0*o55) o57=-XI+XL o58=o57**2 o59=-YI+YL o60=o59**2 o61=-ZI+ZL o62=o61**2 o63=o58+o60+o62 o64=1/o63 o65=o31*o57 o66=o33*o59 o67=o35*o61 o68=o65+o66+o67 o69=o68**2 o70=o22*o38*o49*o52*o56*o64*o69 o71=5.d-1*o5+5.d-1*o6+1.d0*o7 o72=1/o71 o73=o71**2 o74=1/o73 o75=o12*o74 o76=1/(1.d0+1.d0*o75) o77=o30*o43*o49*o64*o69*o72*o76 o78=(1.d0+o2)**2 o79=1/o78 o80=TC**2 o81=1/o80 o82=As*o79*o81 o83=An*(1.d0+As)*o79*o81 o84=An*As*(1.d0+As)*o79*o81 o85=2.5d-1*o82+1.25d-1*o83+1.25d-1*o84 o86=5.d-1*o5+5.d-1*o6+2.5d-1*o7 o87=o86**2 o88=-3.33333333333333d-1*o85+1.111111111111111d-1*o87 o89=sqrt(o88) o90=3.33333333333333d-1*o86+o89 o91=1/o90 o92=o90**2 o93=1/o92 o94=o12*o93 o95=1/(1.d0+2.777777777777778d-2*o94) o96=o21**2 o97=1/o96 o98=o37**4 o99=o49*o97*o98 o100=o29**2 o101=1/o100 o102=o42**4 o103=o101*o102*o49 o104=o63**2 o105=1/o104 o106=o68**4 o107=o105*o106*o49 o108=o103+o107+o99 o109=1/o89 o110=5.d-1*o5-3.33333333333333d-1*o86 o111=o22*o38*o49*o64*o69 o112=3.d0*o103+6.d0*o111 o113=o109*o110*(-1.d0+o112) o114=2.5d-1*o7-3.33333333333333d-1*o86 o115=o30*o43*o49*o64*o69 o116=6.d0*o115+3.d0*o99 o117=o109*o114*(-1.d0+o116) o118=5.d-1*o6-3.33333333333333d-1*o86 o119=o22*o30*o38*o43*o49 o120=3.d0*o107+6.d0*o119 o121=o109*o118*(-1.d0+o120) o122=o113+o117+o121 o123=5.d-1*(-1.d0+3.d0*o108)-1.666666666666667d-1*o122 o124=o123*o91*o95 o125=3.33333333333333d-1*o86-o89 o126=1/o125 o127=o125**2 o128=1/o127 o129=o12*o128 o130=1/(1.d0+2.777777777777778d-2*o129) o131=5.d-1*(-1.d0+3.d0*o108)+1.666666666666667d-1*o122 o132=o126*o130*o131 o133=1.666666666666667d-1*o124+1.666666666666667d-1*o132+6.d0*o & 50+6.d0*o70+6.d0*o77 o134=cc*o133 o135=wh-wn o136=o135**2 o137=o11*o136 o138=1/(1.d0+1.d0*o137) o139=o138*o22*o30*o38*o43*o49*o9 o140=o136*o54 o141=1/(1.d0+1.d0*o140) o142=o141*o22*o38*o49*o52*o64*o69 o143=o136*o74 o144=1/(1.d0+1.d0*o143) o145=o144*o30*o43*o49*o64*o69*o72 o146=o136*o93 o147=1/(1.d0+2.777777777777778d-2*o146) o148=o123*o147*o91 o149=o128*o136 o150=1/(1.d0+2.777777777777778d-2*o149) o151=o126*o131*o150 o152=6.d0*o139+6.d0*o142+6.d0*o145+1.666666666666667d-1*o148+1. & 666666666666667d-1*o151 o153=wh+wn o154=o153**2 o155=o11*o154 o156=1/(1.d0+1.d0*o155) o157=o156*o22*o30*o38*o43*o49*o9 o158=o154*o54 o159=1/(1.d0+1.d0*o158) o160=o159*o22*o38*o49*o52*o64*o69 o161=o154*o74 o162=1/(1.d0+1.d0*o161) o163=o162*o30*o43*o49*o64*o69*o72 o164=o154*o93 o165=1/(1.d0+2.777777777777778d-2*o164) o166=o123*o165*o91 o167=o128*o154 o168=1/(1.d0+2.777777777777778d-2*o167) o169=o126*o131*o168 o170=6.d0*o157+6.d0*o160+6.d0*o163+1.666666666666667d-1*o166+1. & 666666666666667d-1*o169 o171=1.5d0*o133+5.d-1*o152+3.d0*o170 o172=dd*o171 o173=5.d-1*o134+o172 o174=XM-XN o175=o14*o174*o22*o30*o38*o42*o49*o9 o176=o14*o174*o22*o30*o37*o43*o49*o9 o177=o101*o14*o22*o23*o38*o43*o49*o9 o178=o14*o15*o30*o38*o43*o49*o9*o97 o179=o174*o22*o38*o49*o52*o56*o64*o68 o180=o174*o30*o43*o49*o64*o68*o72*o76 o181=o174*o22*o37*o49*o52*o56*o64*o69 o182=o105*o22*o38*o49*o52*o56*o57*o69 o183=o15*o38*o49*o52*o56*o64*o69*o97 o184=o174*o30*o42*o49*o64*o69*o72*o76 o185=o105*o30*o43*o49*o57*o69*o72*o76 o186=o101*o23*o43*o49*o64*o69*o72*o76 o187=o37**3 o188=o174*o187*o49*o97 o189=o21**3 o190=1/o189 o191=o15*o190*o49*o98 o192=o42**3 o193=o101*o174*o192*o49 o194=o29**3 o195=1/o194 o196=o102*o195*o23*o49 o197=o68**3 o198=o105*o174*o197*o49 o199=o63**3 o200=1/o199 o201=o106*o200*o49*o57 o202=4.d0*o188+4.d0*o191+4.d0*o193+4.d0*o196+4.d0*o198+4.d0*o20 & 1 o203=o174*o22*o38*o49*o64*o68 o204=o174*o22*o37*o49*o64*o69 o205=o105*o22*o38*o49*o57*o69 o206=o15*o38*o49*o64*o69*o97 o207=1.2d1*o193+1.2d1*o196+1.2d1*o203+1.2d1*o204+1.2d1*o205+1.2 & d1*o206 o208=o109*o110*o207 o209=o174*o30*o43*o49*o64*o68 o210=o174*o30*o42*o49*o64*o69 o211=o105*o30*o43*o49*o57*o69 o212=o101*o23*o43*o49*o64*o69 o213=1.2d1*o188+1.2d1*o191+1.2d1*o209+1.2d1*o210+1.2d1*o211+1.2 & d1*o212 o214=o109*o114*o213 o215=o174*o22*o30*o38*o42*o49 o216=o174*o22*o30*o37*o43*o49 o217=o101*o22*o23*o38*o43*o49 o218=o15*o30*o38*o43*o49*o97 o219=1.2d1*o198+1.2d1*o201+1.2d1*o215+1.2d1*o216+1.2d1*o217+1.2 & d1*o218 o220=o109*o118*o219 o221=o208+o214+o220 o222=1.5d0*o202-1.666666666666667d-1*o221 o223=o222*o91*o95 o224=1.5d0*o202+1.666666666666667d-1*o221 o225=o126*o130*o224 o226=1.2d1*o175+1.2d1*o176+1.2d1*o177+1.2d1*o178+1.2d1*o179+1.2 & d1*o180+1.2d1*o181+1.2d1*o182+1.2d1*o183+1.2d1*o184+1.2d1*o185+ & 1.2d1*o186+1.666666666666667d-1*o223+1.666666666666667d-1*o225 o227=o138*o174*o22*o30*o38*o42*o49*o9 o228=o138*o174*o22*o30*o37*o43*o49*o9 o229=o101*o138*o22*o23*o38*o43*o49*o9 o230=o138*o15*o30*o38*o43*o49*o9*o97 o231=o141*o174*o22*o38*o49*o52*o64*o68 o232=o144*o174*o30*o43*o49*o64*o68*o72 o233=o141*o174*o22*o37*o49*o52*o64*o69 o234=o105*o141*o22*o38*o49*o52*o57*o69 o235=o141*o15*o38*o49*o52*o64*o69*o97 o236=o144*o174*o30*o42*o49*o64*o69*o72 o237=o105*o144*o30*o43*o49*o57*o69*o72 o238=o101*o144*o23*o43*o49*o64*o69*o72 o239=o147*o222*o91 o240=o126*o150*o224 o241=1.2d1*o227+1.2d1*o228+1.2d1*o229+1.2d1*o230+1.2d1*o231+1.2 & d1*o232+1.2d1*o233+1.2d1*o234+1.2d1*o235+1.2d1*o236+1.2d1*o237+ & 1.2d1*o238+1.666666666666667d-1*o239+1.666666666666667d-1*o240 o242=o156*o174*o22*o30*o38*o42*o49*o9 o243=o156*o174*o22*o30*o37*o43*o49*o9 o244=o101*o156*o22*o23*o38*o43*o49*o9 o245=o15*o156*o30*o38*o43*o49*o9*o97 o246=o159*o174*o22*o38*o49*o52*o64*o68 o247=o162*o174*o30*o43*o49*o64*o68*o72 o248=o159*o174*o22*o37*o49*o52*o64*o69 o249=o105*o159*o22*o38*o49*o52*o57*o69 o250=o15*o159*o38*o49*o52*o64*o69*o97 o251=o162*o174*o30*o42*o49*o64*o69*o72 o252=o105*o162*o30*o43*o49*o57*o69*o72 o253=o101*o162*o23*o43*o49*o64*o69*o72 o254=o165*o222*o91 o255=o126*o168*o224 o256=1.2d1*o242+1.2d1*o243+1.2d1*o244+1.2d1*o245+1.2d1*o246+1.2 & d1*o247+1.2d1*o248+1.2d1*o249+1.2d1*o250+1.2d1*o251+1.2d1*o252+ & 1.2d1*o253+1.666666666666667d-1*o254+1.666666666666667d-1*o255 o257=o173**2 o258=1/o257 o259=YM-YN o260=o14*o22*o259*o30*o38*o42*o49*o9 o261=o14*o22*o259*o30*o37*o43*o49*o9 o262=o101*o14*o22*o25*o38*o43*o49*o9 o263=o14*o17*o30*o38*o43*o49*o9*o97 o264=o22*o259*o38*o49*o52*o56*o64*o68 o265=o259*o30*o43*o49*o64*o68*o72*o76 o266=o22*o259*o37*o49*o52*o56*o64*o69 o267=o105*o22*o38*o49*o52*o56*o59*o69 o268=o17*o38*o49*o52*o56*o64*o69*o97 o269=o259*o30*o42*o49*o64*o69*o72*o76 o270=o105*o30*o43*o49*o59*o69*o72*o76 o271=o101*o25*o43*o49*o64*o69*o72*o76 o272=o187*o259*o49*o97 o273=o17*o190*o49*o98 o274=o101*o192*o259*o49 o275=o102*o195*o25*o49 o276=o105*o197*o259*o49 o277=o106*o200*o49*o59 o278=4.d0*o272+4.d0*o273+4.d0*o274+4.d0*o275+4.d0*o276+4.d0*o27 & 7 o279=o22*o259*o38*o49*o64*o68 o280=o22*o259*o37*o49*o64*o69 o281=o105*o22*o38*o49*o59*o69 o282=o17*o38*o49*o64*o69*o97 o283=1.2d1*o274+1.2d1*o275+1.2d1*o279+1.2d1*o280+1.2d1*o281+1.2 & d1*o282 o284=o109*o110*o283 o285=o259*o30*o43*o49*o64*o68 o286=o259*o30*o42*o49*o64*o69 o287=o105*o30*o43*o49*o59*o69 o288=o101*o25*o43*o49*o64*o69 o289=1.2d1*o272+1.2d1*o273+1.2d1*o285+1.2d1*o286+1.2d1*o287+1.2 & d1*o288 o290=o109*o114*o289 o291=o22*o259*o30*o38*o42*o49 o292=o22*o259*o30*o37*o43*o49 o293=o101*o22*o25*o38*o43*o49 o294=o17*o30*o38*o43*o49*o97 o295=1.2d1*o276+1.2d1*o277+1.2d1*o291+1.2d1*o292+1.2d1*o293+1.2 & d1*o294 o296=o109*o118*o295 o297=o284+o290+o296 o298=1.5d0*o278-1.666666666666667d-1*o297 o299=o298*o91*o95 o300=1.5d0*o278+1.666666666666667d-1*o297 o301=o126*o130*o300 o302=1.2d1*o260+1.2d1*o261+1.2d1*o262+1.2d1*o263+1.2d1*o264+1.2 & d1*o265+1.2d1*o266+1.2d1*o267+1.2d1*o268+1.2d1*o269+1.2d1*o270+ & 1.2d1*o271+1.666666666666667d-1*o299+1.666666666666667d-1*o301 o303=o138*o22*o259*o30*o38*o42*o49*o9 o304=o138*o22*o259*o30*o37*o43*o49*o9 o305=o101*o138*o22*o25*o38*o43*o49*o9 o306=o138*o17*o30*o38*o43*o49*o9*o97 o307=o141*o22*o259*o38*o49*o52*o64*o68 o308=o144*o259*o30*o43*o49*o64*o68*o72 o309=o141*o22*o259*o37*o49*o52*o64*o69 o310=o105*o141*o22*o38*o49*o52*o59*o69 o311=o141*o17*o38*o49*o52*o64*o69*o97 o312=o144*o259*o30*o42*o49*o64*o69*o72 o313=o105*o144*o30*o43*o49*o59*o69*o72 o314=o101*o144*o25*o43*o49*o64*o69*o72 o315=o147*o298*o91 o316=o126*o150*o300 o317=1.2d1*o303+1.2d1*o304+1.2d1*o305+1.2d1*o306+1.2d1*o307+1.2 & d1*o308+1.2d1*o309+1.2d1*o310+1.2d1*o311+1.2d1*o312+1.2d1*o313+ & 1.2d1*o314+1.666666666666667d-1*o315+1.666666666666667d-1*o316 o318=o156*o22*o259*o30*o38*o42*o49*o9 o319=o156*o22*o259*o30*o37*o43*o49*o9 o320=o101*o156*o22*o25*o38*o43*o49*o9 o321=o156*o17*o30*o38*o43*o49*o9*o97 o322=o159*o22*o259*o38*o49*o52*o64*o68 o323=o162*o259*o30*o43*o49*o64*o68*o72 o324=o159*o22*o259*o37*o49*o52*o64*o69 o325=o105*o159*o22*o38*o49*o52*o59*o69 o326=o159*o17*o38*o49*o52*o64*o69*o97 o327=o162*o259*o30*o42*o49*o64*o69*o72 o328=o105*o162*o30*o43*o49*o59*o69*o72 o329=o101*o162*o25*o43*o49*o64*o69*o72 o330=o165*o298*o91 o331=o126*o168*o300 o332=1.2d1*o318+1.2d1*o319+1.2d1*o320+1.2d1*o321+1.2d1*o322+1.2 & d1*o323+1.2d1*o324+1.2d1*o325+1.2d1*o326+1.2d1*o327+1.2d1*o328+ & 1.2d1*o329+1.666666666666667d-1*o330+1.666666666666667d-1*o331 o333=ZM-ZN o334=o14*o22*o30*o333*o38*o42*o49*o9 o335=o14*o22*o30*o333*o37*o43*o49*o9 o336=o101*o14*o22*o27*o38*o43*o49*o9 o337=o14*o19*o30*o38*o43*o49*o9*o97 o338=o22*o333*o38*o49*o52*o56*o64*o68 o339=o30*o333*o43*o49*o64*o68*o72*o76 o340=o22*o333*o37*o49*o52*o56*o64*o69 o341=o105*o22*o38*o49*o52*o56*o61*o69 o342=o19*o38*o49*o52*o56*o64*o69*o97 o343=o30*o333*o42*o49*o64*o69*o72*o76 o344=o105*o30*o43*o49*o61*o69*o72*o76 o345=o101*o27*o43*o49*o64*o69*o72*o76 o346=o187*o333*o49*o97 o347=o19*o190*o49*o98 o348=o101*o192*o333*o49 o349=o102*o195*o27*o49 o350=o105*o197*o333*o49 o351=o106*o200*o49*o61 o352=4.d0*o346+4.d0*o347+4.d0*o348+4.d0*o349+4.d0*o350+4.d0*o35 & 1 o353=o22*o333*o38*o49*o64*o68 o354=o22*o333*o37*o49*o64*o69 o355=o105*o22*o38*o49*o61*o69 o356=o19*o38*o49*o64*o69*o97 o357=1.2d1*o348+1.2d1*o349+1.2d1*o353+1.2d1*o354+1.2d1*o355+1.2 & d1*o356 o358=o109*o110*o357 o359=o30*o333*o43*o49*o64*o68 o360=o30*o333*o42*o49*o64*o69 o361=o105*o30*o43*o49*o61*o69 o362=o101*o27*o43*o49*o64*o69 o363=1.2d1*o346+1.2d1*o347+1.2d1*o359+1.2d1*o360+1.2d1*o361+1.2 & d1*o362 o364=o109*o114*o363 o365=o22*o30*o333*o38*o42*o49 o366=o22*o30*o333*o37*o43*o49 o367=o101*o22*o27*o38*o43*o49 o368=o19*o30*o38*o43*o49*o97 o369=1.2d1*o350+1.2d1*o351+1.2d1*o365+1.2d1*o366+1.2d1*o367+1.2 & d1*o368 o370=o109*o118*o369 o371=o358+o364+o370 o372=1.5d0*o352-1.666666666666667d-1*o371 o373=o372*o91*o95 o374=1.5d0*o352+1.666666666666667d-1*o371 o375=o126*o130*o374 o376=1.2d1*o334+1.2d1*o335+1.2d1*o336+1.2d1*o337+1.2d1*o338+1.2 & d1*o339+1.2d1*o340+1.2d1*o341+1.2d1*o342+1.2d1*o343+1.2d1*o344+ & 1.2d1*o345+1.666666666666667d-1*o373+1.666666666666667d-1*o375 o377=o138*o22*o30*o333*o38*o42*o49*o9 o378=o138*o22*o30*o333*o37*o43*o49*o9 o379=o101*o138*o22*o27*o38*o43*o49*o9 o380=o138*o19*o30*o38*o43*o49*o9*o97 o381=o141*o22*o333*o38*o49*o52*o64*o68 o382=o144*o30*o333*o43*o49*o64*o68*o72 o383=o141*o22*o333*o37*o49*o52*o64*o69 o384=o105*o141*o22*o38*o49*o52*o61*o69 o385=o141*o19*o38*o49*o52*o64*o69*o97 o386=o144*o30*o333*o42*o49*o64*o69*o72 o387=o105*o144*o30*o43*o49*o61*o69*o72 o388=o101*o144*o27*o43*o49*o64*o69*o72 o389=o147*o372*o91 o390=o126*o150*o374 o391=1.2d1*o377+1.2d1*o378+1.2d1*o379+1.2d1*o380+1.2d1*o381+1.2 & d1*o382+1.2d1*o383+1.2d1*o384+1.2d1*o385+1.2d1*o386+1.2d1*o387+ & 1.2d1*o388+1.666666666666667d-1*o389+1.666666666666667d-1*o390 o392=o156*o22*o30*o333*o38*o42*o49*o9 o393=o156*o22*o30*o333*o37*o43*o49*o9 o394=o101*o156*o22*o27*o38*o43*o49*o9 o395=o156*o19*o30*o38*o43*o49*o9*o97 o396=o159*o22*o333*o38*o49*o52*o64*o68 o397=o162*o30*o333*o43*o49*o64*o68*o72 o398=o159*o22*o333*o37*o49*o52*o64*o69 o399=o105*o159*o22*o38*o49*o52*o61*o69 o400=o159*o19*o38*o49*o52*o64*o69*o97 o401=o162*o30*o333*o42*o49*o64*o69*o72 o402=o105*o162*o30*o43*o49*o61*o69*o72 o403=o101*o162*o27*o43*o49*o64*o69*o72 o404=o165*o372*o91 o405=o126*o168*o374 o406=1.2d1*o392+1.2d1*o393+1.2d1*o394+1.2d1*o395+1.2d1*o396+1.2 & d1*o397+1.2d1*o398+1.2d1*o399+1.2d1*o400+1.2d1*o401+1.2d1*o402+ & 1.2d1*o403+1.666666666666667d-1*o404+1.666666666666667d-1*o405 o407=o14*o22*o30*o31*o37*o43*o49*o9 o408=o22*o31*o37*o49*o52*o56*o64*o69 o409=o187*o31*o49*o97 o410=-4.d0*o191+4.d0*o409 o411=-1.2d1*o191+1.2d1*o409 o412=o109*o114*o411 o413=o22*o30*o31*o37*o43*o49 o414=-1.2d1*o218+1.2d1*o413 o415=o109*o118*o414 o416=o22*o31*o37*o49*o64*o69 o417=-1.2d1*o206+1.2d1*o416 o418=o109*o110*o417 o419=o412+o415+o418 o420=1.5d0*o410-1.666666666666667d-1*o419 o421=o420*o91*o95 o422=1.5d0*o410+1.666666666666667d-1*o419 o423=o126*o130*o422 o424=-1.2d1*o178-1.2d1*o183+1.2d1*o407+1.2d1*o408+1.66666666666 & 6667d-1*o421+1.666666666666667d-1*o423 o425=o138*o22*o30*o31*o37*o43*o49*o9 o426=o141*o22*o31*o37*o49*o52*o64*o69 o427=o147*o420*o91 o428=o126*o150*o422 o429=-1.2d1*o230-1.2d1*o235+1.2d1*o425+1.2d1*o426+1.66666666666 & 6667d-1*o427+1.666666666666667d-1*o428 o430=o156*o22*o30*o31*o37*o43*o49*o9 o431=o159*o22*o31*o37*o49*o52*o64*o69 o432=o165*o420*o91 o433=o126*o168*o422 o434=-1.2d1*o245-1.2d1*o250+1.2d1*o430+1.2d1*o431+1.66666666666 & 6667d-1*o432+1.666666666666667d-1*o433 o435=o14*o22*o30*o33*o37*o43*o49*o9 o436=o22*o33*o37*o49*o52*o56*o64*o69 o437=o187*o33*o49*o97 o438=-4.d0*o273+4.d0*o437 o439=-1.2d1*o273+1.2d1*o437 o440=o109*o114*o439 o441=o22*o30*o33*o37*o43*o49 o442=-1.2d1*o294+1.2d1*o441 o443=o109*o118*o442 o444=o22*o33*o37*o49*o64*o69 o445=-1.2d1*o282+1.2d1*o444 o446=o109*o110*o445 o447=o440+o443+o446 o448=1.5d0*o438-1.666666666666667d-1*o447 o449=o448*o91*o95 o450=1.5d0*o438+1.666666666666667d-1*o447 o451=o126*o130*o450 o452=-1.2d1*o263-1.2d1*o268+1.2d1*o435+1.2d1*o436+1.66666666666 & 6667d-1*o449+1.666666666666667d-1*o451 o453=o138*o22*o30*o33*o37*o43*o49*o9 o454=o141*o22*o33*o37*o49*o52*o64*o69 o455=o147*o448*o91 o456=o126*o150*o450 o457=-1.2d1*o306-1.2d1*o311+1.2d1*o453+1.2d1*o454+1.66666666666 & 6667d-1*o455+1.666666666666667d-1*o456 o458=o156*o22*o30*o33*o37*o43*o49*o9 o459=o159*o22*o33*o37*o49*o52*o64*o69 o460=o165*o448*o91 o461=o126*o168*o450 o462=-1.2d1*o321-1.2d1*o326+1.2d1*o458+1.2d1*o459+1.66666666666 & 6667d-1*o460+1.666666666666667d-1*o461 o463=o14*o22*o30*o35*o37*o43*o49*o9 o464=o22*o35*o37*o49*o52*o56*o64*o69 o465=o187*o35*o49*o97 o466=-4.d0*o347+4.d0*o465 o467=-1.2d1*o347+1.2d1*o465 o468=o109*o114*o467 o469=o22*o30*o35*o37*o43*o49 o470=-1.2d1*o368+1.2d1*o469 o471=o109*o118*o470 o472=o22*o35*o37*o49*o64*o69 o473=-1.2d1*o356+1.2d1*o472 o474=o109*o110*o473 o475=o468+o471+o474 o476=1.5d0*o466-1.666666666666667d-1*o475 o477=o476*o91*o95 o478=1.5d0*o466+1.666666666666667d-1*o475 o479=o126*o130*o478 o480=-1.2d1*o337-1.2d1*o342+1.2d1*o463+1.2d1*o464+1.66666666666 & 6667d-1*o477+1.666666666666667d-1*o479 o481=o138*o22*o30*o35*o37*o43*o49*o9 o482=o141*o22*o35*o37*o49*o52*o64*o69 o483=o147*o476*o91 o484=o126*o150*o478 o485=-1.2d1*o380-1.2d1*o385+1.2d1*o481+1.2d1*o482+1.66666666666 & 6667d-1*o483+1.666666666666667d-1*o484 o486=o156*o22*o30*o35*o37*o43*o49*o9 o487=o159*o22*o35*o37*o49*o52*o64*o69 o488=o165*o476*o91 o489=o126*o168*o478 o490=-1.2d1*o395-1.2d1*o400+1.2d1*o486+1.2d1*o487+1.66666666666 & 6667d-1*o488+1.666666666666667d-1*o489 o491=o14*o22*o30*o31*o38*o42*o49*o9 o492=o30*o31*o42*o49*o64*o69*o72*o76 o493=o101*o192*o31*o49 o494=-4.d0*o196+4.d0*o493 o495=o22*o30*o31*o38*o42*o49 o496=-1.2d1*o217+1.2d1*o495 o497=o109*o118*o496 o498=-1.2d1*o196+1.2d1*o493 o499=o109*o110*o498 o500=o30*o31*o42*o49*o64*o69 o501=-1.2d1*o212+1.2d1*o500 o502=o109*o114*o501 o503=o497+o499+o502 o504=1.5d0*o494-1.666666666666667d-1*o503 o505=o504*o91*o95 o506=1.5d0*o494+1.666666666666667d-1*o503 o507=o126*o130*o506 o508=-1.2d1*o177-1.2d1*o186+1.2d1*o491+1.2d1*o492+1.66666666666 & 6667d-1*o505+1.666666666666667d-1*o507 o509=o138*o22*o30*o31*o38*o42*o49*o9 o510=o144*o30*o31*o42*o49*o64*o69*o72 o511=o147*o504*o91 o512=o126*o150*o506 o513=-1.2d1*o229-1.2d1*o238+1.2d1*o509+1.2d1*o510+1.66666666666 & 6667d-1*o511+1.666666666666667d-1*o512 o514=o156*o22*o30*o31*o38*o42*o49*o9 o515=o162*o30*o31*o42*o49*o64*o69*o72 o516=o165*o504*o91 o517=o126*o168*o506 o518=-1.2d1*o244-1.2d1*o253+1.2d1*o514+1.2d1*o515+1.66666666666 & 6667d-1*o516+1.666666666666667d-1*o517 o519=o14*o22*o30*o33*o38*o42*o49*o9 o520=o30*o33*o42*o49*o64*o69*o72*o76 o521=o101*o192*o33*o49 o522=-4.d0*o275+4.d0*o521 o523=o22*o30*o33*o38*o42*o49 o524=-1.2d1*o293+1.2d1*o523 o525=o109*o118*o524 o526=-1.2d1*o275+1.2d1*o521 o527=o109*o110*o526 o528=o30*o33*o42*o49*o64*o69 o529=-1.2d1*o288+1.2d1*o528 o530=o109*o114*o529 o531=o525+o527+o530 o532=1.5d0*o522-1.666666666666667d-1*o531 o533=o532*o91*o95 o534=1.5d0*o522+1.666666666666667d-1*o531 o535=o126*o130*o534 o536=-1.2d1*o262-1.2d1*o271+1.2d1*o519+1.2d1*o520+1.66666666666 & 6667d-1*o533+1.666666666666667d-1*o535 o537=o138*o22*o30*o33*o38*o42*o49*o9 o538=o144*o30*o33*o42*o49*o64*o69*o72 o539=o147*o532*o91 o540=o126*o150*o534 o541=-1.2d1*o305-1.2d1*o314+1.2d1*o537+1.2d1*o538+1.66666666666 & 6667d-1*o539+1.666666666666667d-1*o540 o542=o156*o22*o30*o33*o38*o42*o49*o9 o543=o162*o30*o33*o42*o49*o64*o69*o72 o544=o165*o532*o91 o545=o126*o168*o534 o546=-1.2d1*o320-1.2d1*o329+1.2d1*o542+1.2d1*o543+1.66666666666 & 6667d-1*o544+1.666666666666667d-1*o545 o547=o14*o22*o30*o35*o38*o42*o49*o9 o548=o30*o35*o42*o49*o64*o69*o72*o76 o549=o101*o192*o35*o49 o550=-4.d0*o349+4.d0*o549 o551=o22*o30*o35*o38*o42*o49 o552=-1.2d1*o367+1.2d1*o551 o553=o109*o118*o552 o554=-1.2d1*o349+1.2d1*o549 o555=o109*o110*o554 o556=o30*o35*o42*o49*o64*o69 o557=-1.2d1*o362+1.2d1*o556 o558=o109*o114*o557 o559=o553+o555+o558 o560=1.5d0*o550-1.666666666666667d-1*o559 o561=o560*o91*o95 o562=1.5d0*o550+1.666666666666667d-1*o559 o563=o126*o130*o562 o564=-1.2d1*o336-1.2d1*o345+1.2d1*o547+1.2d1*o548+1.66666666666 & 6667d-1*o561+1.666666666666667d-1*o563 o565=o138*o22*o30*o35*o38*o42*o49*o9 o566=o144*o30*o35*o42*o49*o64*o69*o72 o567=o147*o560*o91 o568=o126*o150*o562 o569=-1.2d1*o379-1.2d1*o388+1.2d1*o565+1.2d1*o566+1.66666666666 & 6667d-1*o567+1.666666666666667d-1*o568 o570=o156*o22*o30*o35*o38*o42*o49*o9 o571=o162*o30*o35*o42*o49*o64*o69*o72 o572=o165*o560*o91 o573=o126*o168*o562 o574=-1.2d1*o394-1.2d1*o403+1.2d1*o570+1.2d1*o571+1.66666666666 & 6667d-1*o572+1.666666666666667d-1*o573 o575=o22*o31*o38*o49*o52*o56*o64*o68 o576=o30*o31*o43*o49*o64*o68*o72*o76 o577=o105*o197*o31*o49 o578=-4.d0*o201+4.d0*o577 o579=o22*o31*o38*o49*o64*o68 o580=-1.2d1*o205+1.2d1*o579 o581=o109*o110*o580 o582=o30*o31*o43*o49*o64*o68 o583=-1.2d1*o211+1.2d1*o582 o584=o109*o114*o583 o585=-1.2d1*o201+1.2d1*o577 o586=o109*o118*o585 o587=o581+o584+o586 o588=1.5d0*o578-1.666666666666667d-1*o587 o589=o588*o91*o95 o590=1.5d0*o578+1.666666666666667d-1*o587 o591=o126*o130*o590 o592=-1.2d1*o182-1.2d1*o185+1.2d1*o575+1.2d1*o576+1.66666666666 & 6667d-1*o589+1.666666666666667d-1*o591 o593=o141*o22*o31*o38*o49*o52*o64*o68 o594=o144*o30*o31*o43*o49*o64*o68*o72 o595=o147*o588*o91 o596=o126*o150*o590 o597=-1.2d1*o234-1.2d1*o237+1.2d1*o593+1.2d1*o594+1.66666666666 & 6667d-1*o595+1.666666666666667d-1*o596 o598=o159*o22*o31*o38*o49*o52*o64*o68 o599=o162*o30*o31*o43*o49*o64*o68*o72 o600=o165*o588*o91 o601=o126*o168*o590 o602=-1.2d1*o249-1.2d1*o252+1.2d1*o598+1.2d1*o599+1.66666666666 & 6667d-1*o600+1.666666666666667d-1*o601 o603=o22*o33*o38*o49*o52*o56*o64*o68 o604=o30*o33*o43*o49*o64*o68*o72*o76 o605=o105*o197*o33*o49 o606=-4.d0*o277+4.d0*o605 o607=o22*o33*o38*o49*o64*o68 o608=-1.2d1*o281+1.2d1*o607 o609=o109*o110*o608 o610=o30*o33*o43*o49*o64*o68 o611=-1.2d1*o287+1.2d1*o610 o612=o109*o114*o611 o613=-1.2d1*o277+1.2d1*o605 o614=o109*o118*o613 o615=o609+o612+o614 o616=1.5d0*o606-1.666666666666667d-1*o615 o617=o616*o91*o95 o618=1.5d0*o606+1.666666666666667d-1*o615 o619=o126*o130*o618 o620=-1.2d1*o267-1.2d1*o270+1.2d1*o603+1.2d1*o604+1.66666666666 & 6667d-1*o617+1.666666666666667d-1*o619 o621=o141*o22*o33*o38*o49*o52*o64*o68 o622=o144*o30*o33*o43*o49*o64*o68*o72 o623=o147*o616*o91 o624=o126*o150*o618 o625=-1.2d1*o310-1.2d1*o313+1.2d1*o621+1.2d1*o622+1.66666666666 & 6667d-1*o623+1.666666666666667d-1*o624 o626=o159*o22*o33*o38*o49*o52*o64*o68 o627=o162*o30*o33*o43*o49*o64*o68*o72 o628=o165*o616*o91 o629=o126*o168*o618 o630=-1.2d1*o325-1.2d1*o328+1.2d1*o626+1.2d1*o627+1.66666666666 & 6667d-1*o628+1.666666666666667d-1*o629 o631=o22*o35*o38*o49*o52*o56*o64*o68 o632=o30*o35*o43*o49*o64*o68*o72*o76 o633=o105*o197*o35*o49 o634=-4.d0*o351+4.d0*o633 o635=o22*o35*o38*o49*o64*o68 o636=-1.2d1*o355+1.2d1*o635 o637=o109*o110*o636 o638=o30*o35*o43*o49*o64*o68 o639=-1.2d1*o361+1.2d1*o638 o640=o109*o114*o639 o641=-1.2d1*o351+1.2d1*o633 o642=o109*o118*o641 o643=o637+o640+o642 o644=1.5d0*o634-1.666666666666667d-1*o643 o645=o644*o91*o95 o646=1.5d0*o634+1.666666666666667d-1*o643 o647=o126*o130*o646 o648=-1.2d1*o341-1.2d1*o344+1.2d1*o631+1.2d1*o632+1.66666666666 & 6667d-1*o645+1.666666666666667d-1*o647 o649=o141*o22*o35*o38*o49*o52*o64*o68 o650=o144*o30*o35*o43*o49*o64*o68*o72 o651=o147*o644*o91 o652=o126*o150*o646 o653=-1.2d1*o384-1.2d1*o387+1.2d1*o649+1.2d1*o650+1.66666666666 & 6667d-1*o651+1.666666666666667d-1*o652 o654=o159*o22*o35*o38*o49*o52*o64*o68 o655=o162*o30*o35*o43*o49*o64*o68*o72 o656=o165*o644*o91 o657=o126*o168*o646 o658=-1.2d1*o399-1.2d1*o402+1.2d1*o654+1.2d1*o655+1.66666666666 & 6667d-1*o656+1.666666666666667d-1*o657 o659=o47**3 o660=1/o659 o661=o14*o22*o30*o31*o38*o43*o660*o9 o662=o22*o31*o38*o52*o56*o64*o660*o69 o663=o30*o31*o43*o64*o660*o69*o72*o76 o664=XI-XK o665=o14*o22*o30*o38*o42*o49*o664*o9 o666=XI-XJ o667=o14*o22*o30*o37*o43*o49*o666*o9 o668=XI-XL o669=o22*o38*o49*o52*o56*o64*o668*o68 o670=o30*o43*o49*o64*o668*o68*o72*o76 o671=o22*o37*o49*o52*o56*o64*o666*o69 o672=o30*o42*o49*o64*o664*o69*o72*o76 o673=o31*o660*o97*o98 o674=o101*o102*o31*o660 o675=o105*o106*o31*o660 o676=o187*o49*o666*o97 o677=o101*o192*o49*o664 o678=o105*o197*o49*o668 o679=4.d0*o673+4.d0*o674+4.d0*o675+4.d0*o676+4.d0*o677+4.d0*o67 & 8 o680=o22*o31*o38*o64*o660*o69 o681=o22*o38*o49*o64*o668*o68 o682=o22*o37*o49*o64*o666*o69 o683=1.2d1*o674+1.2d1*o677+2.4d1*o680+1.2d1*o681+1.2d1*o682 o684=o109*o110*o683 o685=o30*o31*o43*o64*o660*o69 o686=o30*o43*o49*o64*o668*o68 o687=o30*o42*o49*o64*o664*o69 o688=1.2d1*o673+1.2d1*o676+2.4d1*o685+1.2d1*o686+1.2d1*o687 o689=o109*o114*o688 o690=o22*o30*o31*o38*o43*o660 o691=o22*o30*o38*o42*o49*o664 o692=o22*o30*o37*o43*o49*o666 o693=1.2d1*o675+1.2d1*o678+2.4d1*o690+1.2d1*o691+1.2d1*o692 o694=o109*o118*o693 o695=o684+o689+o694 o696=1.5d0*o679-1.666666666666667d-1*o695 o697=o696*o91*o95 o698=1.5d0*o679+1.666666666666667d-1*o695 o699=o126*o130*o698 o700=2.4d1*o661+2.4d1*o662+2.4d1*o663+1.2d1*o665+1.2d1*o667+1.2 & d1*o669+1.2d1*o670+1.2d1*o671+1.2d1*o672+1.666666666666667d-1*o & 697+1.666666666666667d-1*o699 o701=o138*o22*o30*o31*o38*o43*o660*o9 o702=o141*o22*o31*o38*o52*o64*o660*o69 o703=o144*o30*o31*o43*o64*o660*o69*o72 o704=o138*o22*o30*o38*o42*o49*o664*o9 o705=o138*o22*o30*o37*o43*o49*o666*o9 o706=o141*o22*o38*o49*o52*o64*o668*o68 o707=o144*o30*o43*o49*o64*o668*o68*o72 o708=o141*o22*o37*o49*o52*o64*o666*o69 o709=o144*o30*o42*o49*o64*o664*o69*o72 o710=o147*o696*o91 o711=o126*o150*o698 o712=2.4d1*o701+2.4d1*o702+2.4d1*o703+1.2d1*o704+1.2d1*o705+1.2 & d1*o706+1.2d1*o707+1.2d1*o708+1.2d1*o709+1.666666666666667d-1*o & 710+1.666666666666667d-1*o711 o713=o156*o22*o30*o31*o38*o43*o660*o9 o714=o159*o22*o31*o38*o52*o64*o660*o69 o715=o162*o30*o31*o43*o64*o660*o69*o72 o716=o156*o22*o30*o38*o42*o49*o664*o9 o717=o156*o22*o30*o37*o43*o49*o666*o9 o718=o159*o22*o38*o49*o52*o64*o668*o68 o719=o162*o30*o43*o49*o64*o668*o68*o72 o720=o159*o22*o37*o49*o52*o64*o666*o69 o721=o162*o30*o42*o49*o64*o664*o69*o72 o722=o165*o696*o91 o723=o126*o168*o698 o724=2.4d1*o713+2.4d1*o714+2.4d1*o715+1.2d1*o716+1.2d1*o717+1.2 & d1*o718+1.2d1*o719+1.2d1*o720+1.2d1*o721+1.666666666666667d-1*o & 722+1.666666666666667d-1*o723 o725=o14*o22*o30*o33*o38*o43*o660*o9 o726=o22*o33*o38*o52*o56*o64*o660*o69 o727=o30*o33*o43*o64*o660*o69*o72*o76 o728=YI-YK o729=o14*o22*o30*o38*o42*o49*o728*o9 o730=YI-YJ o731=o14*o22*o30*o37*o43*o49*o730*o9 o732=YI-YL o733=o22*o38*o49*o52*o56*o64*o68*o732 o734=o30*o43*o49*o64*o68*o72*o732*o76 o735=o22*o37*o49*o52*o56*o64*o69*o730 o736=o30*o42*o49*o64*o69*o72*o728*o76 o737=o33*o660*o97*o98 o738=o101*o102*o33*o660 o739=o105*o106*o33*o660 o740=o187*o49*o730*o97 o741=o101*o192*o49*o728 o742=o105*o197*o49*o732 o743=4.d0*o737+4.d0*o738+4.d0*o739+4.d0*o740+4.d0*o741+4.d0*o74 & 2 o744=o22*o33*o38*o64*o660*o69 o745=o22*o38*o49*o64*o68*o732 o746=o22*o37*o49*o64*o69*o730 o747=1.2d1*o738+1.2d1*o741+2.4d1*o744+1.2d1*o745+1.2d1*o746 o748=o109*o110*o747 o749=o30*o33*o43*o64*o660*o69 o750=o30*o43*o49*o64*o68*o732 o751=o30*o42*o49*o64*o69*o728 o752=1.2d1*o737+1.2d1*o740+2.4d1*o749+1.2d1*o750+1.2d1*o751 o753=o109*o114*o752 o754=o22*o30*o33*o38*o43*o660 o755=o22*o30*o38*o42*o49*o728 o756=o22*o30*o37*o43*o49*o730 o757=1.2d1*o739+1.2d1*o742+2.4d1*o754+1.2d1*o755+1.2d1*o756 o758=o109*o118*o757 o759=o748+o753+o758 o760=1.5d0*o743-1.666666666666667d-1*o759 o761=o760*o91*o95 o762=1.5d0*o743+1.666666666666667d-1*o759 o763=o126*o130*o762 o764=2.4d1*o725+2.4d1*o726+2.4d1*o727+1.2d1*o729+1.2d1*o731+1.2 & d1*o733+1.2d1*o734+1.2d1*o735+1.2d1*o736+1.666666666666667d-1*o & 761+1.666666666666667d-1*o763 o765=o138*o22*o30*o33*o38*o43*o660*o9 o766=o141*o22*o33*o38*o52*o64*o660*o69 o767=o144*o30*o33*o43*o64*o660*o69*o72 o768=o138*o22*o30*o38*o42*o49*o728*o9 o769=o138*o22*o30*o37*o43*o49*o730*o9 o770=o141*o22*o38*o49*o52*o64*o68*o732 o771=o144*o30*o43*o49*o64*o68*o72*o732 o772=o141*o22*o37*o49*o52*o64*o69*o730 o773=o144*o30*o42*o49*o64*o69*o72*o728 o774=o147*o760*o91 o775=o126*o150*o762 o776=2.4d1*o765+2.4d1*o766+2.4d1*o767+1.2d1*o768+1.2d1*o769+1.2 & d1*o770+1.2d1*o771+1.2d1*o772+1.2d1*o773+1.666666666666667d-1*o & 774+1.666666666666667d-1*o775 o777=o156*o22*o30*o33*o38*o43*o660*o9 o778=o159*o22*o33*o38*o52*o64*o660*o69 o779=o162*o30*o33*o43*o64*o660*o69*o72 o780=o156*o22*o30*o38*o42*o49*o728*o9 o781=o156*o22*o30*o37*o43*o49*o730*o9 o782=o159*o22*o38*o49*o52*o64*o68*o732 o783=o162*o30*o43*o49*o64*o68*o72*o732 o784=o159*o22*o37*o49*o52*o64*o69*o730 o785=o162*o30*o42*o49*o64*o69*o72*o728 o786=o165*o760*o91 o787=o126*o168*o762 o788=2.4d1*o777+2.4d1*o778+2.4d1*o779+1.2d1*o780+1.2d1*o781+1.2 & d1*o782+1.2d1*o783+1.2d1*o784+1.2d1*o785+1.666666666666667d-1*o & 786+1.666666666666667d-1*o787 o789=o14*o22*o30*o35*o38*o43*o660*o9 o790=o22*o35*o38*o52*o56*o64*o660*o69 o791=o30*o35*o43*o64*o660*o69*o72*o76 o792=ZI-ZK o793=o14*o22*o30*o38*o42*o49*o792*o9 o794=ZI-ZJ o795=o14*o22*o30*o37*o43*o49*o794*o9 o796=ZI-ZL o797=o22*o38*o49*o52*o56*o64*o68*o796 o798=o30*o43*o49*o64*o68*o72*o76*o796 o799=o22*o37*o49*o52*o56*o64*o69*o794 o800=o30*o42*o49*o64*o69*o72*o76*o792 o801=o35*o660*o97*o98 o802=o101*o102*o35*o660 o803=o105*o106*o35*o660 o804=o187*o49*o794*o97 o805=o101*o192*o49*o792 o806=o105*o197*o49*o796 o807=4.d0*o801+4.d0*o802+4.d0*o803+4.d0*o804+4.d0*o805+4.d0*o80 & 6 o808=o22*o35*o38*o64*o660*o69 o809=o22*o38*o49*o64*o68*o796 o810=o22*o37*o49*o64*o69*o794 o811=1.2d1*o802+1.2d1*o805+2.4d1*o808+1.2d1*o809+1.2d1*o810 o812=o109*o110*o811 o813=o30*o35*o43*o64*o660*o69 o814=o30*o43*o49*o64*o68*o796 o815=o30*o42*o49*o64*o69*o792 o816=1.2d1*o801+1.2d1*o804+2.4d1*o813+1.2d1*o814+1.2d1*o815 o817=o109*o114*o816 o818=o22*o30*o35*o38*o43*o660 o819=o22*o30*o38*o42*o49*o792 o820=o22*o30*o37*o43*o49*o794 o821=1.2d1*o803+1.2d1*o806+2.4d1*o818+1.2d1*o819+1.2d1*o820 o822=o109*o118*o821 o823=o812+o817+o822 o824=1.5d0*o807-1.666666666666667d-1*o823 o825=o824*o91*o95 o826=1.5d0*o807+1.666666666666667d-1*o823 o827=o126*o130*o826 o828=2.4d1*o789+2.4d1*o790+2.4d1*o791+1.2d1*o793+1.2d1*o795+1.2 & d1*o797+1.2d1*o798+1.2d1*o799+1.2d1*o800+1.666666666666667d-1*o & 825+1.666666666666667d-1*o827 o829=o138*o22*o30*o35*o38*o43*o660*o9 o830=o141*o22*o35*o38*o52*o64*o660*o69 o831=o144*o30*o35*o43*o64*o660*o69*o72 o832=o138*o22*o30*o38*o42*o49*o792*o9 o833=o138*o22*o30*o37*o43*o49*o794*o9 o834=o141*o22*o38*o49*o52*o64*o68*o796 o835=o144*o30*o43*o49*o64*o68*o72*o796 o836=o141*o22*o37*o49*o52*o64*o69*o794 o837=o144*o30*o42*o49*o64*o69*o72*o792 o838=o147*o824*o91 o839=o126*o150*o826 o840=2.4d1*o829+2.4d1*o830+2.4d1*o831+1.2d1*o832+1.2d1*o833+1.2 & d1*o834+1.2d1*o835+1.2d1*o836+1.2d1*o837+1.666666666666667d-1*o & 838+1.666666666666667d-1*o839 o841=o156*o22*o30*o35*o38*o43*o660*o9 o842=o159*o22*o35*o38*o52*o64*o660*o69 o843=o162*o30*o35*o43*o64*o660*o69*o72 o844=o156*o22*o30*o38*o42*o49*o792*o9 o845=o156*o22*o30*o37*o43*o49*o794*o9 o846=o159*o22*o38*o49*o52*o64*o68*o796 o847=o162*o30*o43*o49*o64*o68*o72*o796 o848=o159*o22*o37*o49*o52*o64*o69*o794 o849=o162*o30*o42*o49*o64*o69*o72*o792 o850=o165*o824*o91 o851=o126*o168*o826 o852=2.4d1*o841+2.4d1*o842+2.4d1*o843+1.2d1*o844+1.2d1*o845+1.2 & d1*o846+1.2d1*o847+1.2d1*o848+1.2d1*o849+1.666666666666667d-1*o & 850+1.666666666666667d-1*o851 o853=o14*o22*o23*o30*o38*o42*o49*o9 o854=o14*o15*o22*o30*o37*o43*o49*o9 o855=o22*o38*o49*o52*o56*o57*o64*o68 o856=o30*o43*o49*o57*o64*o68*o72*o76 o857=o15*o22*o37*o49*o52*o56*o64*o69 o858=o23*o30*o42*o49*o64*o69*o72*o76 o859=o15*o187*o49*o97 o860=o101*o192*o23*o49 o861=o105*o197*o49*o57 o862=-4.d0*o673-4.d0*o674-4.d0*o675+4.d0*o859+4.d0*o860+4.d0*o8 & 61 o863=o22*o38*o49*o57*o64*o68 o864=o15*o22*o37*o49*o64*o69 o865=-1.2d1*o674-2.4d1*o680+1.2d1*o860+1.2d1*o863+1.2d1*o864 o866=o109*o110*o865 o867=o30*o43*o49*o57*o64*o68 o868=o23*o30*o42*o49*o64*o69 o869=-1.2d1*o673-2.4d1*o685+1.2d1*o859+1.2d1*o867+1.2d1*o868 o870=o109*o114*o869 o871=o22*o23*o30*o38*o42*o49 o872=o15*o22*o30*o37*o43*o49 o873=-1.2d1*o675-2.4d1*o690+1.2d1*o861+1.2d1*o871+1.2d1*o872 o874=o109*o118*o873 o875=o866+o870+o874 o876=1.5d0*o862-1.666666666666667d-1*o875 o877=o876*o91*o95 o878=1.5d0*o862+1.666666666666667d-1*o875 o879=o126*o130*o878 o880=-2.4d1*o661-2.4d1*o662-2.4d1*o663+1.2d1*o853+1.2d1*o854+1. & 2d1*o855+1.2d1*o856+1.2d1*o857+1.2d1*o858+1.666666666666667d-1* & o877+1.666666666666667d-1*o879 o881=o138*o22*o23*o30*o38*o42*o49*o9 o882=o138*o15*o22*o30*o37*o43*o49*o9 o883=o141*o22*o38*o49*o52*o57*o64*o68 o884=o144*o30*o43*o49*o57*o64*o68*o72 o885=o141*o15*o22*o37*o49*o52*o64*o69 o886=o144*o23*o30*o42*o49*o64*o69*o72 o887=o147*o876*o91 o888=o126*o150*o878 o889=-2.4d1*o701-2.4d1*o702-2.4d1*o703+1.2d1*o881+1.2d1*o882+1. & 2d1*o883+1.2d1*o884+1.2d1*o885+1.2d1*o886+1.666666666666667d-1* & o887+1.666666666666667d-1*o888 o890=o156*o22*o23*o30*o38*o42*o49*o9 o891=o15*o156*o22*o30*o37*o43*o49*o9 o892=o159*o22*o38*o49*o52*o57*o64*o68 o893=o162*o30*o43*o49*o57*o64*o68*o72 o894=o15*o159*o22*o37*o49*o52*o64*o69 o895=o162*o23*o30*o42*o49*o64*o69*o72 o896=o165*o876*o91 o897=o126*o168*o878 o898=-2.4d1*o713-2.4d1*o714-2.4d1*o715+1.2d1*o890+1.2d1*o891+1. & 2d1*o892+1.2d1*o893+1.2d1*o894+1.2d1*o895+1.666666666666667d-1* & o896+1.666666666666667d-1*o897 o899=o14*o22*o25*o30*o38*o42*o49*o9 o900=o14*o17*o22*o30*o37*o43*o49*o9 o901=o22*o38*o49*o52*o56*o59*o64*o68 o902=o30*o43*o49*o59*o64*o68*o72*o76 o903=o17*o22*o37*o49*o52*o56*o64*o69 o904=o25*o30*o42*o49*o64*o69*o72*o76 o905=o17*o187*o49*o97 o906=o101*o192*o25*o49 o907=o105*o197*o49*o59 o908=-4.d0*o737-4.d0*o738-4.d0*o739+4.d0*o905+4.d0*o906+4.d0*o9 & 07 o909=o22*o38*o49*o59*o64*o68 o910=o17*o22*o37*o49*o64*o69 o911=-1.2d1*o738-2.4d1*o744+1.2d1*o906+1.2d1*o909+1.2d1*o910 o912=o109*o110*o911 o913=o30*o43*o49*o59*o64*o68 o914=o25*o30*o42*o49*o64*o69 o915=-1.2d1*o737-2.4d1*o749+1.2d1*o905+1.2d1*o913+1.2d1*o914 o916=o109*o114*o915 o917=o22*o25*o30*o38*o42*o49 o918=o17*o22*o30*o37*o43*o49 o919=-1.2d1*o739-2.4d1*o754+1.2d1*o907+1.2d1*o917+1.2d1*o918 o920=o109*o118*o919 o921=o912+o916+o920 o922=1.5d0*o908-1.666666666666667d-1*o921 o923=o91*o922*o95 o924=1.5d0*o908+1.666666666666667d-1*o921 o925=o126*o130*o924 o926=-2.4d1*o725-2.4d1*o726-2.4d1*o727+1.2d1*o899+1.2d1*o900+1. & 2d1*o901+1.2d1*o902+1.2d1*o903+1.2d1*o904+1.666666666666667d-1* & o923+1.666666666666667d-1*o925 o927=o138*o22*o25*o30*o38*o42*o49*o9 o928=o138*o17*o22*o30*o37*o43*o49*o9 o929=o141*o22*o38*o49*o52*o59*o64*o68 o930=o144*o30*o43*o49*o59*o64*o68*o72 o931=o141*o17*o22*o37*o49*o52*o64*o69 o932=o144*o25*o30*o42*o49*o64*o69*o72 o933=o147*o91*o922 o934=o126*o150*o924 o935=-2.4d1*o765-2.4d1*o766-2.4d1*o767+1.2d1*o927+1.2d1*o928+1. & 2d1*o929+1.2d1*o930+1.2d1*o931+1.2d1*o932+1.666666666666667d-1* & o933+1.666666666666667d-1*o934 o936=o156*o22*o25*o30*o38*o42*o49*o9 o937=o156*o17*o22*o30*o37*o43*o49*o9 o938=o159*o22*o38*o49*o52*o59*o64*o68 o939=o162*o30*o43*o49*o59*o64*o68*o72 o940=o159*o17*o22*o37*o49*o52*o64*o69 o941=o162*o25*o30*o42*o49*o64*o69*o72 o942=o165*o91*o922 o943=o126*o168*o924 o944=-2.4d1*o777-2.4d1*o778-2.4d1*o779+1.2d1*o936+1.2d1*o937+1. & 2d1*o938+1.2d1*o939+1.2d1*o940+1.2d1*o941+1.666666666666667d-1* & o942+1.666666666666667d-1*o943 o945=o14*o22*o27*o30*o38*o42*o49*o9 o946=o14*o19*o22*o30*o37*o43*o49*o9 o947=o22*o38*o49*o52*o56*o61*o64*o68 o948=o30*o43*o49*o61*o64*o68*o72*o76 o949=o19*o22*o37*o49*o52*o56*o64*o69 o950=o27*o30*o42*o49*o64*o69*o72*o76 o951=o187*o19*o49*o97 o952=o101*o192*o27*o49 o953=o105*o197*o49*o61 o954=-4.d0*o801-4.d0*o802-4.d0*o803+4.d0*o951+4.d0*o952+4.d0*o9 & 53 o955=o22*o38*o49*o61*o64*o68 o956=o19*o22*o37*o49*o64*o69 o957=-1.2d1*o802-2.4d1*o808+1.2d1*o952+1.2d1*o955+1.2d1*o956 o958=o109*o110*o957 o959=o30*o43*o49*o61*o64*o68 o960=o27*o30*o42*o49*o64*o69 o961=-1.2d1*o801-2.4d1*o813+1.2d1*o951+1.2d1*o959+1.2d1*o960 o962=o109*o114*o961 o963=o22*o27*o30*o38*o42*o49 o964=o19*o22*o30*o37*o43*o49 o965=-1.2d1*o803-2.4d1*o818+1.2d1*o953+1.2d1*o963+1.2d1*o964 o966=o109*o118*o965 o967=o958+o962+o966 o968=1.5d0*o954-1.666666666666667d-1*o967 o969=o91*o95*o968 o970=1.5d0*o954+1.666666666666667d-1*o967 o971=o126*o130*o970 o972=-2.4d1*o789-2.4d1*o790-2.4d1*o791+1.2d1*o945+1.2d1*o946+1. & 2d1*o947+1.2d1*o948+1.2d1*o949+1.2d1*o950+1.666666666666667d-1* & o969+1.666666666666667d-1*o971 o973=o138*o22*o27*o30*o38*o42*o49*o9 o974=o138*o19*o22*o30*o37*o43*o49*o9 o975=o141*o22*o38*o49*o52*o61*o64*o68 o976=o144*o30*o43*o49*o61*o64*o68*o72 o977=o141*o19*o22*o37*o49*o52*o64*o69 o978=o144*o27*o30*o42*o49*o64*o69*o72 o979=o147*o91*o968 o980=o126*o150*o970 o981=-2.4d1*o829-2.4d1*o830-2.4d1*o831+1.2d1*o973+1.2d1*o974+1. & 2d1*o975+1.2d1*o976+1.2d1*o977+1.2d1*o978+1.666666666666667d-1* & o979+1.666666666666667d-1*o980 o982=o156*o22*o27*o30*o38*o42*o49*o9 o983=o156*o19*o22*o30*o37*o43*o49*o9 o984=o159*o22*o38*o49*o52*o61*o64*o68 o985=o162*o30*o43*o49*o61*o64*o68*o72 o986=o159*o19*o22*o37*o49*o52*o64*o69 o987=o162*o27*o30*o42*o49*o64*o69*o72 o988=o165*o91*o968 o989=o126*o168*o970 o990=-2.4d1*o841-2.4d1*o842-2.4d1*o843+1.2d1*o982+1.2d1*o983+1. & 2d1*o984+1.2d1*o985+1.2d1*o986+1.2d1*o987+1.666666666666667d-1* & o988+1.666666666666667d-1*o989 o991=o22*o30*o38*o43*o49*o9 o992=o22*o38*o49*o52*o64*o69 o993=o30*o43*o49*o64*o69*o72 o994=o123*o91 o995=o126*o131 o996=6.d0*o991+6.d0*o992+6.d0*o993+1.666666666666667d-1*o994+1. & 666666666666667d-1*o995 o997=1.5d0*o133+2.d0*o996 o998=cc*o997 o999=wh**2 o1000=o11*o999 o1001=1/(1.d0+1.d0*o1000) o1002=o1001*o22*o30*o38*o43*o49*o9 o1003=o54*o999 o1004=1/(1.d0+1.d0*o1003) o1005=o1004*o22*o38*o49*o52*o64*o69 o1006=o74*o999 o1007=1/(1.d0+1.d0*o1006) o1008=o1007*o30*o43*o49*o64*o69*o72 o1009=o93*o999 o1010=1/(1.d0+2.777777777777778d-2*o1009) o1011=o1010*o123*o91 o1012=o128*o999 o1013=1/(1.d0+2.777777777777778d-2*o1012) o1014=o1013*o126*o131 o1015=6.d0*o1002+6.d0*o1005+6.d0*o1008+1.666666666666667d-1*o10 & 11+1.666666666666667d-1*o1014 o1016=3.d0*o1015+1.5d0*o133+5.d-1*o152+3.d0*o170+2.d0*o996 o1017=dd*o1016 o1018=5.d-1*o1017+1.666666666666667d-1*o998 o1019=o174*o22*o30*o38*o42*o49*o9 o1020=o174*o22*o30*o37*o43*o49*o9 o1021=o101*o22*o23*o38*o43*o49*o9 o1022=o15*o30*o38*o43*o49*o9*o97 o1023=o174*o22*o38*o49*o52*o64*o68 o1024=o174*o30*o43*o49*o64*o68*o72 o1025=o174*o22*o37*o49*o52*o64*o69 o1026=o105*o22*o38*o49*o52*o57*o69 o1027=o15*o38*o49*o52*o64*o69*o97 o1028=o174*o30*o42*o49*o64*o69*o72 o1029=o105*o30*o43*o49*o57*o69*o72 o1030=o101*o23*o43*o49*o64*o69*o72 o1031=o222*o91 o1032=o126*o224 o1033=1.2d1*o1019+1.2d1*o1020+1.2d1*o1021+1.2d1*o1022+1.2d1*o10 & 23+1.2d1*o1024+1.2d1*o1025+1.2d1*o1026+1.2d1*o1027+1.2d1*o1028+ & 1.2d1*o1029+1.2d1*o1030+1.666666666666667d-1*o1031+1.6666666666 & 66667d-1*o1032 o1034=o1001*o101*o22*o23*o38*o43*o49*o9 o1035=o1001*o15*o30*o38*o43*o49*o9*o97 o1036=o1004*o105*o22*o38*o49*o52*o57*o69 o1037=o1004*o15*o38*o49*o52*o64*o69*o97 o1038=o1007*o105*o30*o43*o49*o57*o69*o72 o1039=o1007*o101*o23*o43*o49*o64*o69*o72 o1040=o1018**2 o1041=1/o1040 o1042=o22*o259*o30*o38*o42*o49*o9 o1043=o22*o259*o30*o37*o43*o49*o9 o1044=o101*o22*o25*o38*o43*o49*o9 o1045=o17*o30*o38*o43*o49*o9*o97 o1046=o22*o259*o38*o49*o52*o64*o68 o1047=o259*o30*o43*o49*o64*o68*o72 o1048=o22*o259*o37*o49*o52*o64*o69 o1049=o105*o22*o38*o49*o52*o59*o69 o1050=o17*o38*o49*o52*o64*o69*o97 o1051=o259*o30*o42*o49*o64*o69*o72 o1052=o105*o30*o43*o49*o59*o69*o72 o1053=o101*o25*o43*o49*o64*o69*o72 o1054=o298*o91 o1055=o126*o300 o1056=1.2d1*o1042+1.2d1*o1043+1.2d1*o1044+1.2d1*o1045+1.2d1*o10 & 46+1.2d1*o1047+1.2d1*o1048+1.2d1*o1049+1.2d1*o1050+1.2d1*o1051+ & 1.2d1*o1052+1.2d1*o1053+1.666666666666667d-1*o1054+1.6666666666 & 66667d-1*o1055 o1057=o1001*o101*o22*o25*o38*o43*o49*o9 o1058=o1001*o17*o30*o38*o43*o49*o9*o97 o1059=o1004*o105*o22*o38*o49*o52*o59*o69 o1060=o1004*o17*o38*o49*o52*o64*o69*o97 o1061=o1007*o105*o30*o43*o49*o59*o69*o72 o1062=o1007*o101*o25*o43*o49*o64*o69*o72 o1063=o22*o30*o333*o38*o42*o49*o9 o1064=o22*o30*o333*o37*o43*o49*o9 o1065=o101*o22*o27*o38*o43*o49*o9 o1066=o19*o30*o38*o43*o49*o9*o97 o1067=o22*o333*o38*o49*o52*o64*o68 o1068=o30*o333*o43*o49*o64*o68*o72 o1069=o22*o333*o37*o49*o52*o64*o69 o1070=o105*o22*o38*o49*o52*o61*o69 o1071=o19*o38*o49*o52*o64*o69*o97 o1072=o30*o333*o42*o49*o64*o69*o72 o1073=o105*o30*o43*o49*o61*o69*o72 o1074=o101*o27*o43*o49*o64*o69*o72 o1075=o372*o91 o1076=o126*o374 o1077=1.2d1*o1063+1.2d1*o1064+1.2d1*o1065+1.2d1*o1066+1.2d1*o10 & 67+1.2d1*o1068+1.2d1*o1069+1.2d1*o1070+1.2d1*o1071+1.2d1*o1072+ & 1.2d1*o1073+1.2d1*o1074+1.666666666666667d-1*o1075+1.6666666666 & 66667d-1*o1076 o1078=o1001*o101*o22*o27*o38*o43*o49*o9 o1079=o1001*o19*o30*o38*o43*o49*o9*o97 o1080=o1004*o105*o22*o38*o49*o52*o61*o69 o1081=o1004*o19*o38*o49*o52*o64*o69*o97 o1082=o1007*o105*o30*o43*o49*o61*o69*o72 o1083=o1007*o101*o27*o43*o49*o64*o69*o72 o1084=o22*o30*o31*o37*o43*o49*o9 o1085=o22*o31*o37*o49*o52*o64*o69 o1086=o420*o91 o1087=o126*o422 o1088=-1.2d1*o1022-1.2d1*o1027+1.2d1*o1084+1.2d1*o1085+1.666666 & 666666667d-1*o1086+1.666666666666667d-1*o1087 o1089=o22*o30*o33*o37*o43*o49*o9 o1090=o22*o33*o37*o49*o52*o64*o69 o1091=o448*o91 o1092=o126*o450 o1093=-1.2d1*o1045-1.2d1*o1050+1.2d1*o1089+1.2d1*o1090+1.666666 & 666666667d-1*o1091+1.666666666666667d-1*o1092 o1094=o22*o30*o35*o37*o43*o49*o9 o1095=o22*o35*o37*o49*o52*o64*o69 o1096=o476*o91 o1097=o126*o478 o1098=-1.2d1*o1066-1.2d1*o1071+1.2d1*o1094+1.2d1*o1095+1.666666 & 666666667d-1*o1096+1.666666666666667d-1*o1097 o1099=o22*o30*o31*o38*o42*o49*o9 o1100=o30*o31*o42*o49*o64*o69*o72 o1101=o504*o91 o1102=o126*o506 o1103=-1.2d1*o1021-1.2d1*o1030+1.2d1*o1099+1.2d1*o1100+1.666666 & 666666667d-1*o1101+1.666666666666667d-1*o1102 o1104=o22*o30*o33*o38*o42*o49*o9 o1105=o30*o33*o42*o49*o64*o69*o72 o1106=o532*o91 o1107=o126*o534 o1108=-1.2d1*o1044-1.2d1*o1053+1.2d1*o1104+1.2d1*o1105+1.666666 & 666666667d-1*o1106+1.666666666666667d-1*o1107 o1109=o22*o30*o35*o38*o42*o49*o9 o1110=o30*o35*o42*o49*o64*o69*o72 o1111=o560*o91 o1112=o126*o562 o1113=-1.2d1*o1065-1.2d1*o1074+1.2d1*o1109+1.2d1*o1110+1.666666 & 666666667d-1*o1111+1.666666666666667d-1*o1112 o1114=o22*o31*o38*o49*o52*o64*o68 o1115=o30*o31*o43*o49*o64*o68*o72 o1116=o588*o91 o1117=o126*o590 o1118=-1.2d1*o1026-1.2d1*o1029+1.2d1*o1114+1.2d1*o1115+1.666666 & 666666667d-1*o1116+1.666666666666667d-1*o1117 o1119=o22*o33*o38*o49*o52*o64*o68 o1120=o30*o33*o43*o49*o64*o68*o72 o1121=o616*o91 o1122=o126*o618 o1123=-1.2d1*o1049-1.2d1*o1052+1.2d1*o1119+1.2d1*o1120+1.666666 & 666666667d-1*o1121+1.666666666666667d-1*o1122 o1124=o22*o35*o38*o49*o52*o64*o68 o1125=o30*o35*o43*o49*o64*o68*o72 o1126=o644*o91 o1127=o126*o646 o1128=-1.2d1*o1070-1.2d1*o1073+1.2d1*o1124+1.2d1*o1125+1.666666 & 666666667d-1*o1126+1.666666666666667d-1*o1127 o1129=o22*o30*o31*o38*o43*o660*o9 o1130=o22*o31*o38*o52*o64*o660*o69 o1131=o30*o31*o43*o64*o660*o69*o72 o1132=o22*o30*o38*o42*o49*o664*o9 o1133=o22*o30*o37*o43*o49*o666*o9 o1134=o22*o38*o49*o52*o64*o668*o68 o1135=o30*o43*o49*o64*o668*o68*o72 o1136=o22*o37*o49*o52*o64*o666*o69 o1137=o30*o42*o49*o64*o664*o69*o72 o1138=o696*o91 o1139=o126*o698 o1140=2.4d1*o1129+2.4d1*o1130+2.4d1*o1131+1.2d1*o1132+1.2d1*o11 & 33+1.2d1*o1134+1.2d1*o1135+1.2d1*o1136+1.2d1*o1137+1.6666666666 & 66667d-1*o1138+1.666666666666667d-1*o1139 o1141=o1001*o22*o30*o31*o38*o43*o660*o9 o1142=o1004*o22*o31*o38*o52*o64*o660*o69 o1143=o1007*o30*o31*o43*o64*o660*o69*o72 o1144=o22*o30*o33*o38*o43*o660*o9 o1145=o22*o33*o38*o52*o64*o660*o69 o1146=o30*o33*o43*o64*o660*o69*o72 o1147=o22*o30*o38*o42*o49*o728*o9 o1148=o22*o30*o37*o43*o49*o730*o9 o1149=o22*o38*o49*o52*o64*o68*o732 o1150=o30*o43*o49*o64*o68*o72*o732 o1151=o22*o37*o49*o52*o64*o69*o730 o1152=o30*o42*o49*o64*o69*o72*o728 o1153=o760*o91 o1154=o126*o762 o1155=2.4d1*o1144+2.4d1*o1145+2.4d1*o1146+1.2d1*o1147+1.2d1*o11 & 48+1.2d1*o1149+1.2d1*o1150+1.2d1*o1151+1.2d1*o1152+1.6666666666 & 66667d-1*o1153+1.666666666666667d-1*o1154 o1156=o1001*o22*o30*o33*o38*o43*o660*o9 o1157=o1004*o22*o33*o38*o52*o64*o660*o69 o1158=o1007*o30*o33*o43*o64*o660*o69*o72 o1159=o22*o30*o35*o38*o43*o660*o9 o1160=o22*o35*o38*o52*o64*o660*o69 o1161=o30*o35*o43*o64*o660*o69*o72 o1162=o22*o30*o38*o42*o49*o792*o9 o1163=o22*o30*o37*o43*o49*o794*o9 o1164=o22*o38*o49*o52*o64*o68*o796 o1165=o30*o43*o49*o64*o68*o72*o796 o1166=o22*o37*o49*o52*o64*o69*o794 o1167=o30*o42*o49*o64*o69*o72*o792 o1168=o824*o91 o1169=o126*o826 o1170=2.4d1*o1159+2.4d1*o1160+2.4d1*o1161+1.2d1*o1162+1.2d1*o11 & 63+1.2d1*o1164+1.2d1*o1165+1.2d1*o1166+1.2d1*o1167+1.6666666666 & 66667d-1*o1168+1.666666666666667d-1*o1169 o1171=o1001*o22*o30*o35*o38*o43*o660*o9 o1172=o1004*o22*o35*o38*o52*o64*o660*o69 o1173=o1007*o30*o35*o43*o64*o660*o69*o72 o1174=o22*o23*o30*o38*o42*o49*o9 o1175=o15*o22*o30*o37*o43*o49*o9 o1176=o22*o38*o49*o52*o57*o64*o68 o1177=o30*o43*o49*o57*o64*o68*o72 o1178=o15*o22*o37*o49*o52*o64*o69 o1179=o23*o30*o42*o49*o64*o69*o72 o1180=o876*o91 o1181=o126*o878 o1182=-2.4d1*o1129-2.4d1*o1130-2.4d1*o1131+1.2d1*o1174+1.2d1*o1 & 175+1.2d1*o1176+1.2d1*o1177+1.2d1*o1178+1.2d1*o1179+1.666666666 & 666667d-1*o1180+1.666666666666667d-1*o1181 o1183=o22*o25*o30*o38*o42*o49*o9 o1184=o17*o22*o30*o37*o43*o49*o9 o1185=o22*o38*o49*o52*o59*o64*o68 o1186=o30*o43*o49*o59*o64*o68*o72 o1187=o17*o22*o37*o49*o52*o64*o69 o1188=o25*o30*o42*o49*o64*o69*o72 o1189=o91*o922 o1190=o126*o924 o1191=-2.4d1*o1144-2.4d1*o1145-2.4d1*o1146+1.2d1*o1183+1.2d1*o1 & 184+1.2d1*o1185+1.2d1*o1186+1.2d1*o1187+1.2d1*o1188+1.666666666 & 666667d-1*o1189+1.666666666666667d-1*o1190 o1192=o22*o27*o30*o38*o42*o49*o9 o1193=o19*o22*o30*o37*o43*o49*o9 o1194=o22*o38*o49*o52*o61*o64*o68 o1195=o30*o43*o49*o61*o64*o68*o72 o1196=o19*o22*o37*o49*o52*o64*o69 o1197=o27*o30*o42*o49*o64*o69*o72 o1198=o91*o968 o1199=o126*o970 o1200=-2.4d1*o1159-2.4d1*o1160-2.4d1*o1161+1.2d1*o1192+1.2d1*o1 & 193+1.2d1*o1194+1.2d1*o1195+1.2d1*o1196+1.2d1*o1197+1.666666666 & 666667d-1*o1198+1.666666666666667d-1*o1199 C C T1 = 1.d0/o173 C C DT1XI = -((5.d-1*cc*o226+dd*(1.5d0*o226+5.d-1*o241+3.d0 & *o256))*o258) DT1YI = -(o258*(5.d-1*cc*o302+dd*(1.5d0*o302+5.d-1*o317+ & 3.d0*o332))) DT1ZI = -(o258*(5.d-1*cc*o376+dd*(1.5d0*o376+5.d-1*o391+ & 3.d0*o406))) DT1XJ = -(o258*(5.d-1*cc*o424+dd*(1.5d0*o424+5.d-1*o429+ & 3.d0*o434))) DT1YJ = -(o258*(5.d-1*cc*o452+dd*(1.5d0*o452+5.d-1*o457+ & 3.d0*o462))) DT1ZJ = -(o258*(5.d-1*cc*o480+dd*(1.5d0*o480+5.d-1*o485+ & 3.d0*o490))) DT1XK = -(o258*(5.d-1*cc*o508+dd*(1.5d0*o508+5.d-1*o513+ & 3.d0*o518))) DT1YK = -(o258*(5.d-1*cc*o536+dd*(1.5d0*o536+5.d-1*o541+ & 3.d0*o546))) DT1ZK = -(o258*(5.d-1*cc*o564+dd*(1.5d0*o564+5.d-1*o569+ & 3.d0*o574))) DT1XL = -(o258*(5.d-1*cc*o592+dd*(1.5d0*o592+5.d-1*o597+ & 3.d0*o602))) DT1YL = -(o258*(5.d-1*cc*o620+dd*(1.5d0*o620+5.d-1*o625+ & 3.d0*o630))) DT1ZL = -(o258*(5.d-1*cc*o648+dd*(1.5d0*o648+5.d-1*o653+ & 3.d0*o658))) DT1XM = -(o258*(5.d-1*cc*o700+dd*(1.5d0*o700+5.d-1*o712+ & 3.d0*o724))) DT1YM = -(o258*(5.d-1*cc*o764+dd*(1.5d0*o764+5.d-1*o776+ & 3.d0*o788))) DT1ZM = -(o258*(5.d-1*cc*o828+dd*(1.5d0*o828+5.d-1*o840+ & 3.d0*o852))) DT1XN = -(o258*(5.d-1*cc*o880+dd*(1.5d0*o880+5.d-1*o889+ & 3.d0*o898))) DT1YN = -(o258*(5.d-1*cc*o926+dd*(1.5d0*o926+5.d-1*o935+ & 3.d0*o944))) DT1ZN = -(o258*(5.d-1*cc*o972+dd*(1.5d0*o972+5.d-1*o981+ & 3.d0*o990))) C T2 = 1.d0/o1018 C C DT2XI = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1033+1.5d0* & o226)+5.d-1* & dd*(2.d0*o1033+1.5d0*o226+5.d-1*o241+3.d0*o256+3.d0*(1.2d1*o103 & 4+1.2d1*o1035+1.2d1*o1036+1.2d1*o1037+1.2d1*o1038+1.2d1*o1039+1 & .666666666666667d-1*o1013*o126*o224+1.2d1*o1004*o174*o22*o38*o4 & 9*o52*o64*o68+1.2d1*o1004*o174*o22*o37*o49*o52*o64*o69+1.2d1*o1 & 007*o174*o30*o43*o49*o64*o68*o72+1.2d1*o1007*o174*o30*o42*o49*o & 64*o69*o72+1.2d1*o1001*o174*o22*o30*o38*o42*o49*o9+1.2d1*o1001* & o174*o22*o30*o37*o43*o49*o9+1.666666666666667d-1*o1010*o222*o91 & )))) DT2YI = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1056+1.5d0* & o302)+5.d-1* & dd*(2.d0*o1056+1.5d0*o302+5.d-1*o317+3.d0*o332+3.d0*(1.2d1*o105 & 7+1.2d1*o1058+1.2d1*o1059+1.2d1*o1060+1.2d1*o1061+1.2d1*o1062+1 & .666666666666667d-1*o1013*o126*o300+1.2d1*o1004*o22*o259*o38*o4 & 9*o52*o64*o68+1.2d1*o1004*o22*o259*o37*o49*o52*o64*o69+1.2d1*o1 & 007*o259*o30*o43*o49*o64*o68*o72+1.2d1*o1007*o259*o30*o42*o49*o & 64*o69*o72+1.2d1*o1001*o22*o259*o30*o38*o42*o49*o9+1.2d1*o1001* & o22*o259*o30*o37*o43*o49*o9+1.666666666666667d-1*o1010*o298*o91 & )))) DT2ZI = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1077+1.5d0* & o376)+5.d-1* & dd*(2.d0*o1077+1.5d0*o376+5.d-1*o391+3.d0*o406+3.d0*(1.2d1*o107 & 8+1.2d1*o1079+1.2d1*o1080+1.2d1*o1081+1.2d1*o1082+1.2d1*o1083+1 & .666666666666667d-1*o1013*o126*o374+1.2d1*o1004*o22*o333*o38*o4 & 9*o52*o64*o68+1.2d1*o1004*o22*o333*o37*o49*o52*o64*o69+1.2d1*o1 & 007*o30*o333*o43*o49*o64*o68*o72+1.2d1*o1007*o30*o333*o42*o49*o & 64*o69*o72+1.2d1*o1001*o22*o30*o333*o38*o42*o49*o9+1.2d1*o1001* & o22*o30*o333*o37*o43*o49*o9+1.666666666666667d-1*o1010*o372*o91 & )))) DT2XJ = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1088+1.5d0* & o424)+5.d-1* & dd*(2.d0*o1088+1.5d0*o424+5.d-1*o429+3.d0*o434+3.d0*(-1.2d1*o10 & 35-1.2d1*o1037+1.666666666666667d-1*o1013*o126*o422+1.2d1*o1004 & *o22*o31*o37*o49*o52*o64*o69+1.2d1*o1001*o22*o30*o31*o37*o43*o4 & 9*o9+1.666666666666667d-1*o1010*o420*o91)))) DT2YJ = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1093+1.5d0* & o452)+5.d-1* & dd*(2.d0*o1093+1.5d0*o452+5.d-1*o457+3.d0*o462+3.d0*(-1.2d1*o10 & 58-1.2d1*o1060+1.666666666666667d-1*o1013*o126*o450+1.2d1*o1004 & *o22*o33*o37*o49*o52*o64*o69+1.2d1*o1001*o22*o30*o33*o37*o43*o4 & 9*o9+1.666666666666667d-1*o1010*o448*o91)))) DT2ZJ = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1098+1.5d0* & o480)+5.d-1* & dd*(2.d0*o1098+1.5d0*o480+5.d-1*o485+3.d0*o490+3.d0*(-1.2d1*o10 & 79-1.2d1*o1081+1.666666666666667d-1*o1013*o126*o478+1.2d1*o1004 & *o22*o35*o37*o49*o52*o64*o69+1.2d1*o1001*o22*o30*o35*o37*o43*o4 & 9*o9+1.666666666666667d-1*o1010*o476*o91)))) DT2XK = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1103+1.5d0* & o508)+5.d-1* & dd*(2.d0*o1103+1.5d0*o508+5.d-1*o513+3.d0*o518+3.d0*(-1.2d1*o10 & 34-1.2d1*o1039+1.666666666666667d-1*o1013*o126*o506+1.2d1*o1007 & *o30*o31*o42*o49*o64*o69*o72+1.2d1*o1001*o22*o30*o31*o38*o42*o4 & 9*o9+1.666666666666667d-1*o1010*o504*o91)))) DT2YK = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1108+1.5d0* & o536)+5.d-1* & dd*(2.d0*o1108+1.5d0*o536+5.d-1*o541+3.d0*o546+3.d0*(-1.2d1*o10 & 57-1.2d1*o1062+1.666666666666667d-1*o1013*o126*o534+1.2d1*o1007 & *o30*o33*o42*o49*o64*o69*o72+1.2d1*o1001*o22*o30*o33*o38*o42*o4 & 9*o9+1.666666666666667d-1*o1010*o532*o91)))) DT2ZK = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1113+1.5d0* & o564)+5.d-1* & dd*(2.d0*o1113+1.5d0*o564+5.d-1*o569+3.d0*o574+3.d0*(-1.2d1*o10 & 78-1.2d1*o1083+1.666666666666667d-1*o1013*o126*o562+1.2d1*o1007 & *o30*o35*o42*o49*o64*o69*o72+1.2d1*o1001*o22*o30*o35*o38*o42*o4 & 9*o9+1.666666666666667d-1*o1010*o560*o91)))) DT2XL = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1118+1.5d0* & o592)+5.d-1* & dd*(2.d0*o1118+1.5d0*o592+5.d-1*o597+3.d0*o602+3.d0*(-1.2d1*o10 & 36-1.2d1*o1038+1.666666666666667d-1*o1013*o126*o590+1.2d1*o1004 & *o22*o31*o38*o49*o52*o64*o68+1.2d1*o1007*o30*o31*o43*o49*o64*o6 & 8*o72+1.666666666666667d-1*o1010*o588*o91)))) DT2YL = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1123+1.5d0* & o620)+5.d-1* & dd*(2.d0*o1123+1.5d0*o620+5.d-1*o625+3.d0*o630+3.d0*(-1.2d1*o10 & 59-1.2d1*o1061+1.666666666666667d-1*o1013*o126*o618+1.2d1*o1004 & *o22*o33*o38*o49*o52*o64*o68+1.2d1*o1007*o30*o33*o43*o49*o64*o6 & 8*o72+1.666666666666667d-1*o1010*o616*o91)))) DT2ZL = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1128+1.5d0* & o648)+5.d-1* & dd*(2.d0*o1128+1.5d0*o648+5.d-1*o653+3.d0*o658+3.d0*(-1.2d1*o10 & 80-1.2d1*o1082+1.666666666666667d-1*o1013*o126*o646+1.2d1*o1004 & *o22*o35*o38*o49*o52*o64*o68+1.2d1*o1007*o30*o35*o43*o49*o64*o6 & 8*o72+1.666666666666667d-1*o1010*o644*o91)))) DT2XM = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1140+1.5d0* & o700)+5.d-1* & dd*(2.d0*o1140+1.5d0*o700+5.d-1*o712+3.d0*o724+3.d0*(2.4d1*o114 & 1+2.4d1*o1142+2.4d1*o1143+1.2d1*o1004*o22*o38*o49*o52*o64*o668* & o68+1.2d1*o1004*o22*o37*o49*o52*o64*o666*o69+1.666666666666667d & -1*o1013*o126*o698+1.2d1*o1007*o30*o43*o49*o64*o668*o68*o72+1.2 & d1*o1007*o30*o42*o49*o64*o664*o69*o72+1.2d1*o1001*o22*o30*o38*o & 42*o49*o664*o9+1.2d1*o1001*o22*o30*o37*o43*o49*o666*o9+1.666666 & 666666667d-1*o1010*o696*o91)))) DT2YM = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1155+1.5d0* & o764)+5.d-1* & dd*(2.d0*o1155+1.5d0*o764+5.d-1*o776+3.d0*o788+3.d0*(2.4d1*o115 & 6+2.4d1*o1157+2.4d1*o1158+1.2d1*o1007*o30*o42*o49*o64*o69*o72*o & 728+1.2d1*o1004*o22*o37*o49*o52*o64*o69*o730+1.2d1*o1004*o22*o3 & 8*o49*o52*o64*o68*o732+1.2d1*o1007*o30*o43*o49*o64*o68*o72*o732 & +1.666666666666667d-1*o1013*o126*o762+1.2d1*o1001*o22*o30*o38*o & 42*o49*o728*o9+1.2d1*o1001*o22*o30*o37*o43*o49*o730*o9+1.666666 & 666666667d-1*o1010*o760*o91)))) DT2ZM = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1170+1.5d0* & o828)+5.d-1* & dd*(2.d0*o1170+1.5d0*o828+5.d-1*o840+3.d0*o852+3.d0*(2.4d1*o117 & 1+2.4d1*o1172+2.4d1*o1173+1.2d1*o1007*o30*o42*o49*o64*o69*o72*o & 792+1.2d1*o1004*o22*o37*o49*o52*o64*o69*o794+1.2d1*o1004*o22*o3 & 8*o49*o52*o64*o68*o796+1.2d1*o1007*o30*o43*o49*o64*o68*o72*o796 & +1.666666666666667d-1*o1013*o126*o826+1.2d1*o1001*o22*o30*o38*o & 42*o49*o792*o9+1.2d1*o1001*o22*o30*o37*o43*o49*o794*o9+1.666666 & 666666667d-1*o1010*o824*o91)))) DT2XN = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1182+1.5d0* & o880)+5.d-1* & dd*(2.d0*o1182+1.5d0*o880+5.d-1*o889+3.d0*o898+3.d0*(-2.4d1*o11 & 41-2.4d1*o1142-2.4d1*o1143+1.2d1*o1004*o22*o38*o49*o52*o57*o64* & o68+1.2d1*o1004*o15*o22*o37*o49*o52*o64*o69+1.2d1*o1007*o30*o43 & *o49*o57*o64*o68*o72+1.2d1*o1007*o23*o30*o42*o49*o64*o69*o72+1. & 666666666666667d-1*o1013*o126*o878+1.2d1*o1001*o22*o23*o30*o38* & o42*o49*o9+1.2d1*o1001*o15*o22*o30*o37*o43*o49*o9+1.66666666666 & 6667d-1*o1010*o876*o91)))) DT2YN = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1191+1.5d0* & o926)+5.d-1* & dd*(2.d0*o1191+3.d0*(-2.4d1*o1156-2.4d1*o1157-2.4d1*o1158+1.2d1 & *o1004*o22*o38*o49*o52*o59*o64*o68+1.2d1*o1004*o17*o22*o37*o49* & o52*o64*o69+1.2d1*o1007*o30*o43*o49*o59*o64*o68*o72+1.2d1*o1007 & *o25*o30*o42*o49*o64*o69*o72+1.2d1*o1001*o22*o25*o30*o38*o42*o4 & 9*o9+1.2d1*o1001*o17*o22*o30*o37*o43*o49*o9+1.666666666666667d- & 1*o1010*o91*o922+1.666666666666667d-1*o1013*o126*o924)+1.5d0*o9 & 26+5.d-1*o935+3.d0*o944))) DT2ZN = -(o1041*(1.666666666666667d-1*cc*(2.d0*o1200+1.5d0* & o972)+5.d-1* & dd*(2.d0*o1200+3.d0*(-2.4d1*o1171-2.4d1*o1172-2.4d1*o1173+1.2d1 & *o1004*o22*o38*o49*o52*o61*o64*o68+1.2d1*o1004*o19*o22*o37*o49* & o52*o64*o69+1.2d1*o1007*o30*o43*o49*o61*o64*o68*o72+1.2d1*o1007 & *o27*o30*o42*o49*o64*o69*o72+1.2d1*o1001*o22*o27*o30*o38*o42*o4 & 9*o9+1.2d1*o1001*o19*o22*o30*o37*o43*o49*o9+1.666666666666667d- & 1*o1010*o91*o968+1.666666666666667d-1*o1013*o126*o970)+1.5d0*o9 & 72+5.d-1*o981+3.d0*o990))) C C C C C calculate energy as a function of delta anis C C C C Casting the coefficients to relate to the correct definition C COEF0 = Correlation time Tc C COEF1 = Dz / 0.5*(Dx+Dy) C COEF2 = 1.5*(Dy-Dx)/(Dz-0.5*(Dy+Dx)) C COEF3 = Wh C COEF4 = Wn C C C CALCDANI = T1/T2 IF (WHICH.EQ.'ANALYZE') DANICALC(COUNT) = CALCDANI DELTADANI = CALCDANI-OBSDANI IF ( (ABS(DELTADANI).LT.ERRDANI) & .AND.(POTENTIAL.EQ.SQUARE)) THEN E = 0 DF1 = 0 DF2 = 0 ELSE E = K1*(DELTADANI**2) DF1 = 2*K1*DELTADANI*(1.0/T2) DF2 = 2*K1*DELTADANI*(-1.0*T1/(T2**2)) END IF C C C accumulate energy C EDA = EDA + E C C now update forces if in energy & force mode C IF (WHICH.NE.'ANALYZE') THEN DX(ATOMI(COUNT))=DX(ATOMI(COUNT))+DF1*DT1XI+DF2*DT2XI DY(ATOMI(COUNT))=DY(ATOMI(COUNT))+DF1*DT1YI+DF2*DT2YI DZ(ATOMI(COUNT))=DZ(ATOMI(COUNT))+DF1*DT1ZI+DF2*DT2ZI DX(ATOMJ(COUNT))=DX(ATOMJ(COUNT))+DF1*DT1XJ+DF2*DT2XJ DY(ATOMJ(COUNT))=DY(ATOMJ(COUNT))+DF1*DT1YJ+DF2*DT2YJ DZ(ATOMJ(COUNT))=DZ(ATOMJ(COUNT))+DF1*DT1ZJ+DF2*DT2ZJ DX(ATOMK(COUNT))=DX(ATOMK(COUNT))+DF1*DT1XK+DF2*DT2XK DY(ATOMK(COUNT))=DY(ATOMK(COUNT))+DF1*DT1YK+DF2*DT2YK DZ(ATOMK(COUNT))=DZ(ATOMK(COUNT))+DF1*DT1ZK+DF2*DT2ZK DX(ATOML(COUNT))=DX(ATOML(COUNT))+DF1*DT1XL+DF2*DT2XL DY(ATOML(COUNT))=DY(ATOML(COUNT))+DF1*DT1YL+DF2*DT2YL DZ(ATOML(COUNT))=DZ(ATOML(COUNT))+DF1*DT1ZL+DF2*DT2ZL DX(ATOMM(COUNT))=DX(ATOMM(COUNT))+DF1*DT1XM+DF2*DT2XM DY(ATOMM(COUNT))=DY(ATOMM(COUNT))+DF1*DT1YM+DF2*DT2YM DZ(ATOMM(COUNT))=DZ(ATOMM(COUNT))+DF1*DT1ZM+DF2*DT2ZM DX(ATOMN(COUNT))=DX(ATOMN(COUNT))+DF1*DT1XN+DF2*DT2XN DY(ATOMN(COUNT))=DY(ATOMN(COUNT))+DF1*DT1YN+DF2*DT2YN DZ(ATOMN(COUNT))=DZ(ATOMN(COUNT))+DF1*DT1ZN+DF2*DT2ZN END IF END DO RETURN END C C================ SUBROUTINE READDANI C C reads in coupling constant information C C by Nico Tjandra June 1997 C================ IMPLICIT NONE C include files INCLUDE 'cns.inc' INCLUDE 'comand.inc' INCLUDE 'diff_anis.inc' INCLUDE 'funct.inc' INCLUDE 'mtf.inc' INCLUDE 'heap.inc' INCLUDE 'numbers.inc' C i/o C local variables INTEGER COUNT, SPTR, OLDCLASS, OLDMAXDANI DOUBLE PRECISION K1, CUTOFF, COEF0, COEF1, COEF2, COEF3, COEF4 CHARACTER*4 THENAME C begin C C this is used by READDANI2 to hold the selection C SPTR=ALLHP(INTEG4(NATOM)) C C reset database only if no couplings have been entered C ie., this is the first time in the cns script that C DANIsotropy has appeared C IF (DANIIPTR.EQ.0) THEN CALL DANIDEFAULTS CALL ALLOCDANI(0, MAXDANI) END IF C C now read input C CALL PUSEND('DANI>') DO WHILE (.NOT.DONE) CALL NEXTWD('DANI>') CALL MISCOM('DANI>',USED) IF (.NOT.USED) THEN C IF (WD(1:4).EQ.'HELP') THEN C CALL CNSHELP('cns-danisotropy') C C Get class name. Determine if it's an already-defined class. C Insert a new class if it's not. C ELSE IF (WD(1:4).EQ.'CLAS') THEN OLDCLASS = CURCLASS CALL NEXTA4('class name =', THENAME) MODE = NEW DO COUNT = 1, NCLASSES IF (DANICLASSNAMES(COUNT).EQ.THENAME) THEN MODE = UPDATE CURCLASS = COUNT END IF END DO IF (MODE.EQ.NEW) THEN C C make sure you can't add more than the maximum C number of classes C IF (OLDCLASS.EQ.MAXDANICLASSES) THEN CALL DSPERR('DANI','Too many classes.') CALL DSPERR('DANI', & 'Increase MAXDANICLASSES and recompile.') CALL WRNDIE(-5, 'READDANI', & 'Too many Anisotropy classes.') END IF NCLASSES = NCLASSES + 1 CURCLASS = NCLASSES DANICLASSNAMES(CURCLASS) = THENAME C C If this isn't the first class, close off the old class C IF (NCLASSES.GT.1) THEN DANIASSNDX(OLDCLASS) = NDANI END IF C END IF C C C set force constant for current class C ELSE IF (WD(1:4).EQ.'FORC') THEN CALL NEXTF('force constant =', K1) C C start a default class if there isn't one defined C IF (CURCLASS.EQ.0) THEN NCLASSES = 1 CURCLASS = 1 END IF WRITE(PUNIT, '(A, A, A, F8.3)') & 'Setting force const for class ', & DANICLASSNAMES(CURCLASS), ' to ', K1 DANIFORCES(CURCLASS) = K1 C C set coefficient constant for current class C ELSE IF (WD(1:4).EQ.'COEF') THEN CALL NEXTF('coefficient 0 =', COEF0) CALL NEXTF('coefficient 1 =', COEF1) CALL NEXTF('coefficient 2 =', COEF2) CALL NEXTF('coefficient 3 =', COEF3) CALL NEXTF('coefficient 4 =', COEF4) C C start a default class if there isn't one already defined C IF (CURCLASS.EQ.0) THEN NCLASSES = 1 CURCLASS = 1 END IF WRITE(PUNIT, '(A,A,A, F8.3, F8.3, F8.3, F8.3, F8.3)') & 'Setting coefficients for class ', & DANICLASSNAMES(CURCLASS), 'to', COEF0, COEF1, COEF2, & COEF3, COEF4 C DANICOEF0(CURCLASS) = COEF0 DANICOEF1(CURCLASS) = COEF1 DANICOEF2(CURCLASS) = COEF2 DANICOEF3(CURCLASS) = COEF3 DANICOEF4(CURCLASS) = COEF4 C C C reset couplings database C ELSE IF (WD(1:4).EQ.'RESE') THEN CALL DANIHP CALL DANIINIT CALL ALLOCDANI(0, MAXDANI) C C set potential type C ELSE IF (WD(1:4).EQ.'POTE') THEN CALL NEXTA4('potential type =', THENAME) IF (THENAME.EQ.'SQUA') THEN WRITE(PUNIT, '(A)') 'using square well potential.' POTENTIAL = SQUARE ELSE IF (THENAME.EQ.'HARM') THEN WRITE(PUNIT, '(A)') 'using harmonic potential.' POTENTIAL = HARMONIC ELSE CALL DSPERR('DANI','unknown potential. Using square.') POTENTIAL = SQUARE END IF C C C change number of assignment slots C ELSE IF (WD(1:4).EQ.'NRES') THEN OLDMAXDANI = MAXDANI CALL NEXTI('number of slots =', MAXDANI) CALL ALLOCDANI(OLDMAXDANI, MAXDANI) WRITE(PUNIT, '(A, I8)') & 'DANI: Allocating space for', MAXDANI, & 'number of constraints ' C C C read in an assignment C ELSE IF (WD(1:4).EQ.'ASSI') THEN C C make sure you can't add more coupling assignments C than you have slots for C IF (NDANI.EQ.MAXDANI) THEN CALL DSPERR('DANI','Too many assignments.') CALL DSPERR('DANI', & 'Increase NREStraints and run again.') CALL WRNDIE(-1,'DANI>', & 'exceeded allocation for DANI restraints') END IF C C if there isn't a class specified, C start a default class C IF (CURCLASS.EQ.0) THEN NCLASSES = 1 CURCLASS = 1 END IF CALL READDANI2(HEAP(DANIIPTR), HEAP(DANIJPTR), & HEAP(DANIKPTR), HEAP(DANILPTR), & HEAP(DANIMPTR), & HEAP(DANINPTR), HEAP(DANIOBSPTR), & HEAP(DANIERRPTR), HEAP(SPTR)) C C print violations C ELSE IF (WD(1:4).EQ.'PRIN') THEN CALL NEXTWD('PRINt>') IF (WD(1:4).NE.'THRE') THEN CALL DSPERR('DANI', & 'print expects THREshold parameter.') ELSE CALL NEXTF('THREshold =', CUTOFF) IF (CUTOFF.LT.ZERO) THEN CALL DSPERR('DANI', 'cutoff must be positive.') CUTOFF = ABS(CUTOFF) END IF CALL NEXTA4('ALL or CLASs>', THENAME) IF (THENAME(1:3).EQ.'ALL') THEN DO COUNT = 1,NCLASSES PRINTCLASS(COUNT) = .TRUE. END DO ELSE IF (THENAME(1:4).EQ.'CLAS') THEN CALL NEXTA4('class name =', THENAME) DO COUNT = 1,NCLASSES IF (DANICLASSNAMES(COUNT).EQ. & THENAME) THEN PRINTCLASS(COUNT) = .TRUE. ELSE PRINTCLASS(COUNT) = .FALSE. END IF END DO ELSE DO COUNT = 1,NCLASSES PRINTCLASS(COUNT) = .TRUE. END DO END IF C CALL PRINTDANI(CUTOFF, HEAP(CALCDANIPTR), HEAP(DANIOBSPTR), & HEAP(DANIIPTR), HEAP(DANIJPTR), HEAP(DANIKPTR), & HEAP(DANILPTR), HEAP(DANIMPTR), HEAP(DANINPTR)) END IF C C check for END statement C ELSE CALL CHKEND('DANI>', DONE) END IF END IF END DO DONE = .FALSE. CALL FREHP(SPTR,INTEG4(NATOM)) RETURN END C=============== SUBROUTINE ALLOCDANI (OLDSIZE, NEWSIZE) C C resets coupling constant arrays to hold SIZE entries C C by Nico Tjandra June 1997 C=============== IMPLICIT NONE C include files INCLUDE 'funct.inc' INCLUDE 'diff_anis.inc' C i/o INTEGER OLDSIZE, NEWSIZE C begin IF (OLDSIZE.NE.0) THEN CALL FREHP(DANIIPTR, INTEG4(OLDSIZE)) CALL FREHP(DANIJPTR, INTEG4(OLDSIZE)) CALL FREHP(DANIKPTR, INTEG4(OLDSIZE)) CALL FREHP(DANILPTR, INTEG4(OLDSIZE)) CALL FREHP(DANIMPTR, INTEG4(OLDSIZE)) CALL FREHP(DANINPTR, INTEG4(OLDSIZE)) CALL FREHP(DANIOBSPTR, IREAL8(OLDSIZE)) CALL FREHP(DANIERRPTR, IREAL8(OLDSIZE)) CALL FREHP(CALCDANIPTR, IREAL8(OLDSIZE)) END IF C DANIIPTR = 0 DANIJPTR = 0 DANIKPTR = 0 DANILPTR = 0 DANIMPTR = 0 DANINPTR = 0 DANIOBSPTR = 0 DANIERRPTR = 0 CALCDANIPTR = 0 C IF (NEWSIZE.NE.0) THEN DANIIPTR = ALLHP(INTEG4(NEWSIZE)) DANIJPTR = ALLHP(INTEG4(NEWSIZE)) DANIKPTR = ALLHP(INTEG4(NEWSIZE)) DANILPTR = ALLHP(INTEG4(NEWSIZE)) DANIMPTR = ALLHP(INTEG4(NEWSIZE)) DANINPTR = ALLHP(INTEG4(NEWSIZE)) DANIOBSPTR = ALLHP(IREAL8(NEWSIZE)) DANIERRPTR = ALLHP(IREAL8(NEWSIZE)) CALCDANIPTR = ALLHP(IREAL8(NEWSIZE)) END IF C RETURN END C============== SUBROUTINE DANIDEFAULTS C C sets up defaults C C by Nico Tjandra June 1997 C============== IMPLICIT NONE C include files INCLUDE 'diff_anis.inc' C local variables INTEGER COUNT C begin DANIIPTR=0 MODE = NEW MAXDANI = 200 NDANI = 0 NCLASSES = 0 CURCLASS = 0 POTENTIAL = SQUARE DO COUNT = 1, MAXDANICLASSES DANICLASSNAMES(COUNT) = 'DEFAULT' DANIASSNDX(COUNT) = 0 DANIFORCES (COUNT) = 50 DANICOEF0(COUNT) = 1.0 DANICOEF1(COUNT) = 1.0 DANICOEF2(COUNT) = 1.0 DANICOEF3(COUNT) = 1.0 DANICOEF4(COUNT) = 1.0 END DO RETURN END C================= SUBROUTINE READDANI2 (ATOMI, ATOMJ, ATOMK, ATOML, & ATOMM, ATOMN, DANIOBS, DANIERR, SEL) C C reads actual J coupling assignments into arrays C C by Nico Tjandra June 1997 C============== IMPLICIT NONE C include files INCLUDE 'cns.inc' INCLUDE 'coord.inc' INCLUDE 'diff_anis.inc' INCLUDE 'mtf.inc' C i/o INTEGER ATOMI(*), ATOMJ(*), ATOMK(*), ATOML(*) INTEGER ATOMM(*), ATOMN(*), SEL(*) DOUBLE PRECISION DANIOBS(*), DANIERR(*) C local variables INTEGER NSEL, INSERTPOS, COUNT, CURSTOP, OTHERSTOP, NFLAG DOUBLE PRECISION DANIO, DANIE C begin C C if we're in update mode, make a space for the new line C NFLAG = 0 C IF (MODE.EQ.UPDATE) THEN DO COUNT = NDANI+1, DANIASSNDX(CURCLASS)+1, -1 ATOMI(COUNT) = ATOMI(COUNT-1) ATOMJ(COUNT) = ATOMJ(COUNT-1) ATOMK(COUNT) = ATOMK(COUNT-1) ATOML(COUNT) = ATOML(COUNT-1) ATOMM(COUNT) = ATOMM(COUNT-1) ATOMN(COUNT) = ATOMN(COUNT-1) DANIOBS(COUNT) = DANIOBS(COUNT-1) DANIERR(COUNT) = DANIERR(COUNT-1) END DO CURSTOP = DANIASSNDX(CURCLASS) DO COUNT = 1, NCLASSES OTHERSTOP = DANIASSNDX(COUNT) IF (OTHERSTOP.GT.CURSTOP) THEN DANIASSNDX(COUNT) = OTHERSTOP + 1 END IF END DO DANIASSNDX(CURCLASS) = CURSTOP + 1 INSERTPOS = CURSTOP NDANI = NDANI + 1 ELSE NDANI = NDANI + 1 INSERTPOS = NDANI DANIASSNDX(CURCLASS) = INSERTPOS END IF C CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.) IF (NSEL.GT.1) THEN CALL DSPERR('DANI', & 'more than 1 atom in sel for atom i. Using first') END IF IF (NSEL.EQ.0) THEN CALL DSPERR('DANI','Error with atom selection') NFLAG = 1 ENDIF CALL MAKIND(SEL, NATOM, NSEL) ATOMI(INSERTPOS) = SEL(1) C CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.) IF (NSEL.GT.1) THEN CALL DSPERR('DANI', & 'more than 1 atom in sel for atom j. Using first') END IF IF (NSEL.EQ.0) THEN CALL DSPERR('DANI','Error with atom selection') NFLAG = 1 ENDIF CALL MAKIND(SEL, NATOM, NSEL) ATOMJ(INSERTPOS) = SEL(1) C CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.) IF (NSEL.GT.1) THEN CALL DSPERR('DANI', & 'more than 1 atom in sel for atom k. Using first') END IF IF (NSEL.EQ.0) THEN CALL DSPERR('DANI','Error with atom selection') NFLAG = 1 ENDIF CALL MAKIND(SEL, NATOM, NSEL) ATOMK(INSERTPOS) = SEL(1) C CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.) IF (NSEL.GT.1) THEN CALL DSPERR('DANI', & 'more than 1 atom in sel for atom l. Using first') END IF IF (NSEL.EQ.0) THEN CALL DSPERR('DANI','Error with atom selection') NFLAG = 1 ENDIF CALL MAKIND(SEL, NATOM, NSEL) ATOML(INSERTPOS) = SEL(1) C CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.) IF (NSEL.GT.1) THEN CALL DSPERR('DANI', & 'more than 1 atom in sel for atom l. Using first') END IF IF (NSEL.EQ.0) THEN CALL DSPERR('DANI','Error with atom selection') NFLAG = 1 ENDIF CALL MAKIND(SEL, NATOM, NSEL) ATOMM(INSERTPOS) = SEL(1) C CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.) IF (NSEL.GT.1) THEN CALL DSPERR('DANI', & 'more than 1 atom in sel for atom l. Using first') END IF IF (NSEL.EQ.0) THEN CALL DSPERR('DANI','Error with atom selection') NFLAG = 1 ENDIF CALL MAKIND(SEL, NATOM, NSEL) ATOMN(INSERTPOS) = SEL(1) C CALL NEXTF('observed Anis =', DANIO) CALL NEXTF('error in Anis =', DANIE) DANIOBS(INSERTPOS) =DANIO DANIERR(INSERTPOS) = DANIE C C Check for error atom selection. If there is one than reset the C counter for restraint C IF (NFLAG.EQ.1) THEN NDANI = NDANI-1 ENDIF C RETURN END C=============== SUBROUTINE DANIINIT C C initializes J-coupling stuff C C by nico Tjandra June 1997 C================ IMPLICIT NONE C include files INCLUDE 'diff_anis.inc' C begin CALL DANIDEFAULTS RETURN END C=============== SUBROUTINE DANIHP C C deallocates diffusion anisotropy stuff C C by Gregory Warren 1/9/98 C================ IMPLICIT NONE C include files INCLUDE 'diff_anis.inc' C begin IF (DANIIPTR.NE.0) THEN CALL ALLOCDANI(MAXDANI,0) END IF RETURN END C================= SUBROUTINE PRINTDANI (CUTOFF, DANICALC, DANIOBS, & ATOMI, ATOMJ, ATOMK, ATOML, ATOMM, ATOMN) C C prints anisotropy with deltaJ greater than cutoff C calculates RMS deviation and puts it into $RESULT C C by Nico Tjandra June 1997 C================= IMPLICIT NONE C include files INCLUDE 'cns.inc' INCLUDE 'diff_anis.inc' INCLUDE 'comand.inc' INCLUDE 'mtf.inc' INCLUDE 'numbers.inc' C i/o DOUBLE PRECISION CUTOFF, DANICALC(*), DANIOBS(*) INTEGER ATOMI(*), ATOMJ(*), ATOMK(*), ATOML(*) INTEGER ATOMM(*), ATOMN(*) C local variables DOUBLE PRECISION CALCDANI, OBSDANI, DELTADANI, DANIENERGY INTEGER COUNT, CLASS, I, J, K, L, M, N DOUBLE PRECISION RMS, VIOLS, DBPREC DOUBLE COMPLEX DUMMY2 LOGICAL PRINTTHISCLASS C begin RMS = ZERO VIOLS = ZERO C C make sure that the calcJ array is up to date C CALL EDANI(DANIENERGY, 'ANALYZE') WRITE (PUNIT, '(A)') 'The following anisotropies have' WRITE (PUNIT, '(A)') 'delta DANI greater than the cutoff:' WRITE (PUNIT, '(A)') & '(calc Diff-Anis) (obs Diff-Anis) (delta Diff-Anis)' C C write out first class heading C CLASS = 1 PRINTTHISCLASS = PRINTCLASS(CLASS) IF (PRINTTHISCLASS) THEN WRITE (PUNIT, '(A, A)') 'class ', DANICLASSNAMES(1) END IF C C for every anisotropy entry, C COUNT = 0 C DO WHILE(COUNT.LT.NDANI) COUNT=COUNT+1 C C is this the start of a new class? C IF (DANIASSNDX(CLASS).LT.COUNT) THEN CLASS = CLASS + 1 PRINTTHISCLASS = PRINTCLASS(CLASS) IF (DANIASSNDX(CLASS).EQ.DANIASSNDX(CLASS-1)) THEN COUNT = COUNT-1 END IF IF (PRINTTHISCLASS) THEN WRITE (PUNIT, '(A, A)') 'class ', DANICLASSNAMES(CLASS) END IF END IF C C If this assignment is in a class to be printed IF((PRINTTHISCLASS).AND.(DANIASSNDX(CLASS).NE. & DANIASSNDX(CLASS-1))) THEN C C always update RMS delta J C CALCDANI = DANICALC(COUNT) OBSDANI = DANIOBS(COUNT) DELTADANI = (CALCDANI-OBSDANI) RMS = RMS + DELTADANI**2 C C print out delta Js greater than cutoff C and update number of violations C IF (ABS(DELTADANI).GT.CUTOFF) THEN I = ATOMI(COUNT) J = ATOMJ(COUNT) K = ATOMK(COUNT) L = ATOML(COUNT) M = ATOMM(COUNT) N = ATOMN(COUNT) WRITE(PUNIT,'(5X,8(1X,A), 3(F8.3))') & SEGID(M),RESID(M),RES(M),TYPE(M), & SEGID(N),RESID(N),RES(N),TYPE(N), & CALCDANI, OBSDANI, DELTADANI VIOLS = VIOLS + ONE END IF END IF END DO C IF (NDANI.GT.0) THEN RMS = SQRT(RMS / NDANI) ELSE RMS = 0.0 ENDIF DBPREC = NDANI CALL DECLAR('NUMBER', 'DP', ' ', DUMMY2, DBPREC) CALL DECLAR('RMS', 'DP', ' ', DUMMY2, RMS) CALL DECLAR('VIOLATIONS', 'DP', ' ', DUMMY2, VIOLS) RETURN END C SUBROUTINE SCRADANI C C Author: Gregory Warren IMPLICIT NONE C I/O INCLUDE 'diff_anis.inc' C reset diffusion anisotropy database IF (DANIIPTR.NE.0) THEN WRITE(6,'(A)') & ' SCRATC-warning: diffusion anistropy database erased.' CALL DANIHP CALL DANIINIT END IF RETURN END C