(* Kimball model and energy derived from DNA Facio PDB ES 10 April 2013 *) 
Clear[tw,n2,n,nb,bc,bh,u,u1,u2,Vee,Vne,Vnn];

t1=1.889726*{{-0.808,-8.873,-2.080},{-0.196,-10.198,-1.860},{-2.029,-8.583,-1.290},{0.270,-7.725,-1.830},{4.189,-7.682,-0.250},{2.252,-5.410,-0.470},{3.108,-4.880,-0.052},{2.035,-6.740,0.240},{2.420,-6.732,1.260},{0.982,-7.015,0.289},{2.818,-7.701,-0.640},{2.437,-8.720,-0.571},{2.608,-7.128,-2.040},{3.470,-7.426,-2.637},{1.364,-7.578,-2.770},{1.104,-6.836,-3.525},{1.561,-8.532,-3.258},{2.530,-5.682,-1.830},{-1.645,-0.802,-0.140},{-1.755,0.217,-0.042},{-2.476,-1.404,-0.230},{-0.417,-1.347,-0.150},{-0.191,-2.733,-0.280},{-1.049,-3.808,-0.430},{-0.253,-4.833,-0.500},{-0.608,-5.858,-0.612},{1.073,-4.504,-0.420},{1.112,-3.139,-0.280},{2.210,-2.370,-0.160},{1.868,-1.105,-0.040},{2.698,-0.407,0.067},{0.658,-0.554,-0.030},{4.562,-7.653,1.300},{5.836,-8.366,1.520},{3.403,-8.136,2.090},{4.759,-6.091,1.550},{7.904,-3.753,3.130},{5.002,-3.053,2.910},{5.383,-2.121,3.328},{5.608,-4.257,3.620},{5.915,-4.024,4.640},{4.918,-5.099,3.669},{6.806,-4.574,2.740},{7.097,-5.622,2.809},{6.300,-4.234,1.340},{7.172,-3.968,0.743},{5.558,-5.329,0.610},{4.911,-4.882,-0.145},{6.278,-5.985,0.122},{5.387,-3.110,1.550},{-0.508,-2.712,3.110},{-1.100,-3.551,3.036},{-0.940,-1.785,3.223},{0.810,-2.827,3.060},{1.548,-1.714,3.160},{2.902,-1.771,3.110},{3.611,-0.761,3.200},{3.515,-3.013,2.960},{2.775,-4.147,2.860},{3.281,-5.104,2.737},{1.416,-4.113,2.910},{0.819,-5.022,2.839},{8.189,-3.510,4.680},{9.638,-3.338,4.900},{7.536,-4.582,5.470},{7.430,-2.130,4.930},{8.600,1.610,6.510},{5.841,0.470,6.290},{5.602,1.448,6.708},{7.039,-0.147,7.000},{7.150,0.221,8.020},{6.975,-1.234,7.049},{8.195,0.300,6.120},{9.046,-0.377,6.189},{7.585,0.278,4.720},{8.134,1.006,4.123},{7.629,-1.044,3.990},{6.843,-1.063,3.235},{8.597,-1.152,3.502},{6.186,0.650,4.930},{0.254,-1.812,6.620},{-0.749,-1.602,6.718},{0.570,-2.789,6.530},{1.152,-0.813,6.610},{2.540,-1.026,6.480},{3.297,-2.174,6.330},{4.518,-1.734,6.260},{5.383,-2.388,6.148},{4.615,-0.371,6.340},{3.329,0.088,6.480},{2.937,1.369,6.600},{1.628,1.435,6.720},{1.221,2.440,6.827},{0.730,0.455,6.730},{8.688,1.973,8.060},{9.759,2.965,8.280},{8.790,0.723,8.850},{7.263,2.644,8.310},{6.012,6.358,9.890},{4.449,3.814,9.670},{3.681,4.464,10.088},{5.781,4.018,10.380},{5.655,4.382,11.400},{6.369,3.101,10.429},{6.453,5.060,9.500},{7.540,5.013,9.569},{5.973,4.683,8.100},{5.990,5.595,7.503},{6.786,3.639,7.370},{6.160,3.162,6.615},{7.632,4.121,6.882},{4.622,4.162,8.310},{1.271,-1.317,10.000},{0.336,-1.736,10.098},{2.100,-1.921,9.910},{1.410,0.020,9.990},{2.658,0.663,9.860},{3.946,0.179,9.710},{4.675,1.253,9.640},{5.759,1.232,9.528},{3.952,2.412,9.720},{2.642,2.028,9.860},{1.571,2.834,9.980},{0.474,2.118,10.100},{-0.447,2.692,10.207},{0.324,0.797,10.110},{5.869,6.703,11.440},{6.153,8.135,11.660},{6.686,5.751,12.230},{4.322,6.408,11.690},{1.126,8.677,13.270},{1.358,5.700,13.050},{0.354,5.775,13.468},{2.315,6.649,13.760},{1.999,6.869,14.780},{3.329,6.252,13.809},{2.247,7.887,12.880},{3.154,8.487,12.949},{2.080,7.300,11.480},{1.557,8.047,10.883},{3.351,6.932,10.750},{3.125,6.180,9.995},{3.752,7.820,10.262},{1.293,6.084,11.690},{-2.299,1.943,13.620},{-2.894,2.784,13.627},{-2.726,1.010,13.713},{1.681,-0.315,13.390},{1.143,0.792,13.370},{1.735,2.068,13.240},{3.077,2.430,13.100},{3.046,3.735,13.020},{3.943,4.343,12.901},{1.779,4.275,13.100},{0.934,3.196,13.240},{-0.418,3.263,13.360},{-0.974,2.061,13.490},{-0.255,0.884,13.490},{-0.783,0.005,13.585},{0.808,8.873,14.820},{0.196,10.198,15.040},{2.029,8.583,15.610},{-0.270,7.725,15.070},{-4.189,7.682,16.650},{-2.252,5.410,16.430},{-3.108,4.880,16.848},{-2.035,6.740,17.140},{-2.420,6.732,18.160},{-0.982,7.015,17.189},{-2.818,7.701,16.260},{-2.437,8.720,16.329},{-2.608,7.128,14.860},{-3.470,7.426,14.263},{-1.364,7.578,14.130},{-1.104,6.836,13.375},{-1.561,8.532,13.642},{-2.530,5.682,15.070},{1.645,0.802,16.760},{1.755,-0.217,16.858},{2.476,1.404,16.670},{0.417,1.347,16.750},{0.191,2.733,16.620},{1.049,3.808,16.470},{0.253,4.833,16.400},{0.608,5.858,16.288},{-1.073,4.504,16.480},{-1.112,3.139,16.620},{-2.210,2.370,16.740},{-1.868,1.105,16.860},{-2.698,0.407,16.967},{-0.658,0.554,16.870},{-4.562,7.653,18.200},{-5.836,8.366,18.420},{-3.403,8.136,18.990},{-4.759,6.091,18.450},{-7.904,3.753,20.030},{-5.002,3.053,19.810},{-5.383,2.121,20.228},{-5.608,4.257,20.520},{-5.915,4.024,21.540},{-4.918,5.099,20.569},{-6.806,4.574,19.640},{-7.097,5.622,19.709},{-6.300,4.234,18.240},{-7.172,3.968,17.643},{-5.558,5.329,17.510},{-4.911,4.882,16.755},{-6.278,5.985,17.022},{-5.387,3.110,18.450},{0.859,1.616,20.140},{1.548,0.856,20.238},{1.178,2.591,20.050},{-0.454,1.335,20.130},{-1.452,2.323,20.000},{-1.390,3.697,19.850},{-2.636,4.059,19.780},{-2.951,5.096,19.668},{-3.515,3.013,19.860},{-2.745,1.886,20.000},{-3.181,0.618,20.120},{-2.161,-0.204,20.240},{-2.422,-1.257,20.347},{-0.858,0.061,20.250},{-8.189,3.510,21.580},{-9.638,3.338,21.800},{-7.536,4.582,22.370},{-7.430,2.130,21.830},{-8.600,-1.610,23.410},{-5.841,-0.470,23.190},{-5.602,-1.448,23.608},{-7.039,0.147,23.900},{-7.150,-0.221,24.920},{-6.975,1.234,23.949},{-8.195,-0.300,23.020},{-9.046,0.377,23.089},{-7.585,-0.278,21.620},{-8.134,-1.006,21.023},{-7.629,1.044,20.890},{-6.843,1.063,20.135},{-8.597,1.152,20.402},{-6.186,-0.650,21.830},{-1.184,2.493,23.390},{-1.197,3.520,23.316},{-0.288,1.996,23.503},{-2.317,1.811,23.340},{-2.260,0.477,23.440},{-3.389,-0.273,23.390},{-3.369,-1.507,23.480},{-4.615,0.371,23.240},{-4.683,1.723,23.140},{-5.655,2.201,23.017},{-3.563,2.495,23.190},{-3.614,3.582,23.119},{-8.688,-1.973,24.960},{-9.759,-2.965,25.180},{-8.790,-0.723,25.750},{-7.263,-2.644,25.210},{-6.012,-6.358,26.790},{-4.449,-3.814,26.570},{-3.681,-4.464,26.988},{-5.781,-4.018,27.280},{-5.655,-4.382,28.300},{-6.369,-3.101,27.329},{-6.453,-5.060,26.400},{-7.540,-5.013,26.469},{-5.973,-4.683,25.000},{-5.990,-5.595,24.403},{-6.786,-3.639,24.270},{-6.160,-3.162,23.515},{-7.632,-4.121,23.782},{-4.622,-4.162,25.210},{-2.483,1.337,26.770},{-2.973,0.203,26.720},{-2.176,-0.914,26.810},{-1.165,-0.753,26.918},{-2.593,-2.230,26.770},{-1.809,-3.159,26.860},{-3.952,-2.412,26.620},{-4.823,-1.356,26.520},{-5.889,-1.548,26.398},{-4.379,-0.084,26.570},{-5.286,1.105,26.460},{-4.689,2.015,26.420},{-5.884,1.023,25.552},{-5.946,1.141,27.327},{-5.869,-6.703,28.340},{-6.153,-8.135,28.560},{-6.686,-5.751,29.130},{-4.322,-6.408,28.590},{-1.126,-8.677,30.170},{-1.358,-5.700,29.950},{-0.354,-5.775,30.368},{-2.315,-6.649,30.660},{-1.999,-6.869,31.680},{-3.329,-6.252,30.709},{-2.247,-7.887,29.780},{-3.154,-8.487,29.849},{-2.080,-7.300,28.380},{-1.557,-8.047,27.783},{-3.351,-6.932,27.650},{-3.125,-6.180,26.895},{-3.752,-7.820,27.162},{-1.293,-6.084,28.590},{-1.802,0.318,30.280},{-1.292,1.208,30.378},{-2.828,0.320,30.190},{-1.129,-0.845,30.270},{-1.761,-2.099,30.140},{-3.087,-2.464,29.990},{-3.045,-3.761,29.920},{-3.935,-4.382,29.808},{-1.779,-4.275,30.000},{-0.945,-3.193,30.140},{0.395,-3.216,30.260},{0.862,-1.992,30.380},{1.944,-1.915,30.487},{0.207,-0.835,30.390},{4.562,7.653,32.500},{5.836,8.366,32.280},{3.403,8.136,31.710},{4.759,6.091,32.250},{7.904,3.753,30.670},{5.002,3.053,30.890},{5.383,2.121,30.472},{5.608,4.257,30.180},{5.915,4.024,29.160},{4.918,5.099,30.131},{6.806,4.574,31.060},{7.097,5.622,30.991},{6.300,4.234,32.460},{7.172,3.968,33.057},{5.558,5.329,33.190},{4.911,4.882,33.945},{6.278,5.985,33.679},{5.387,3.110,32.250},{-0.505,2.774,30.690},{0.726,2.891,30.740},{1.542,1.787,30.650},{1.076,0.875,30.542},{2.922,1.777,30.690},{3.563,0.745,30.600},{3.515,3.013,30.840},{2.780,4.168,30.940},{3.292,5.122,31.062},{1.434,4.139,30.890},{0.583,5.369,31.000},{1.221,6.252,31.014},{-0.092,5.422,30.146},{0.001,5.327,31.920},{8.189,3.510,29.120},{9.638,3.338,28.900},{7.536,4.582,28.330},{7.430,2.130,28.870},{8.600,-1.610,27.290},{5.841,-0.470,27.510},{5.602,-1.448,27.092},{7.039,0.147,26.800},{7.150,-0.221,25.780},{6.975,1.234,26.751},{8.195,-0.300,27.680},{9.046,0.377,27.611},{7.585,-0.278,29.080},{8.134,-1.006,29.677},{7.629,1.045,29.810},{6.843,1.063,30.565},{8.597,1.152,30.299},{6.186,-0.650,28.870},{0.254,1.812,27.180},{-0.749,1.602,27.082},{0.570,2.789,27.270},{1.152,0.813,27.190},{2.540,1.026,27.320},{3.297,2.174,27.470},{4.518,1.734,27.540},{5.383,2.388,27.652},{4.615,0.371,27.460},{3.329,-0.088,27.320},{2.937,-1.369,27.200},{1.628,-1.435,27.080},{1.221,-2.440,26.973},{0.730,-0.455,27.070},{8.688,-1.973,25.740},{9.759,-2.965,25.520},{8.790,-0.723,24.950},{7.263,-2.644,25.490},{6.012,-6.358,23.910},{4.449,-3.814,24.130},{3.681,-4.464,23.712},{5.781,-4.018,23.420},{5.655,-4.382,22.400},{6.369,-3.101,23.371},{6.453,-5.060,24.300},{7.540,-5.012,24.231},{5.973,-4.683,25.700},{5.990,-5.595,26.297},{6.786,-3.639,26.430},{6.160,-3.162,27.185},{7.632,-4.121,26.919},{4.622,-4.162,25.490},{-0.718,-2.923,23.560},{-0.705,-3.953,23.553},{-1.612,-2.419,23.467},{1.175,1.243,23.790},{1.390,0.031,23.810},{2.619,-0.653,23.940},{3.917,-0.157,24.080},{4.660,-1.231,24.160},{5.742,-1.196,24.279},{3.952,-2.412,24.080},{2.635,-2.037,23.940},{1.580,-2.885,23.820},{0.423,-2.240,23.690},{0.313,-0.865,23.690},{-0.631,-0.464,23.595},{5.869,-6.703,22.360},{6.153,-8.135,22.140},{6.686,-5.751,21.570},{4.322,-6.408,22.110},{1.126,-8.677,20.530},{1.358,-5.700,20.750},{0.354,-5.775,20.332},{2.315,-6.649,20.040},{1.999,-6.869,19.020},{3.329,-6.252,19.991},{2.247,-7.887,20.920},{3.154,-8.487,20.851},{2.080,-7.300,22.320},{1.557,-8.047,22.917},{3.351,-6.932,23.050},{3.125,-6.180,23.805},{3.752,-7.820,23.539},{1.293,-6.084,22.110},{2.795,-0.377,20.550},{2.525,-1.583,20.600},{1.223,-2.019,20.510},{0.500,-1.293,20.402},{0.787,-3.328,20.550},{-0.393,-3.619,20.460},{1.779,-4.274,20.700},{3.105,-3.932,20.800},{3.854,-4.714,20.922},{3.493,-2.642,20.750},{4.926,-2.213,20.860},{5.569,-3.093,20.874},{5.185,-1.588,20.006},{5.066,-1.647,21.780},{0.808,-8.873,18.980},{0.196,-10.198,18.760},{2.029,-8.583,18.190},{-0.270,-7.725,18.730},{-4.189,-7.682,17.150},{-2.252,-5.410,17.370},{-3.108,-4.880,16.952},{-2.035,-6.740,16.660},{-2.420,-6.732,15.640},{-0.982,-7.015,16.611},{-2.818,-7.701,17.540},{-2.437,-8.720,17.471},{-2.608,-7.128,18.940},{-3.470,-7.426,19.537},{-1.364,-7.578,19.670},{-1.104,-6.836,20.425},{-1.561,-8.532,20.159},{-2.530,-5.682,18.730},{2.039,-1.948,17.170},{1.112,-2.765,17.220},{-0.197,-2.352,17.130},{-0.356,-1.340,17.022},{-1.320,-3.155,17.170},{-2.445,-2.697,17.080},{-1.073,-4.504,17.320},{0.201,-5.006,17.420},{0.347,-6.079,17.542},{1.273,-4.191,17.370},{2.684,-4.686,17.480},{2.687,-5.776,17.494},{3.261,-4.332,16.626},{3.131,-4.310,18.400},{-4.562,-7.653,15.600},{-5.836,-8.366,15.380},{-3.403,-8.136,14.810},{-4.759,-6.091,15.350},{-7.904,-3.753,13.770},{-5.002,-3.053,13.990},{-5.383,-2.121,13.572},{-5.608,-4.257,13.280},{-5.915,-4.024,12.260},{-4.918,-5.099,13.231},{-6.806,-4.574,14.160},{-7.097,-5.622,14.091},{-6.300,-4.234,15.560},{-7.172,-3.968,16.157},{-5.558,-5.329,16.290},{-4.911,-4.882,17.045},{-6.278,-5.985,16.779},{-5.387,-3.110,15.350},{0.508,-2.712,13.790},{1.100,-3.551,13.864},{0.940,-1.785,13.677},{-0.810,-2.827,13.840},{-1.548,-1.714,13.740},{-2.902,-1.771,13.790},{-3.611,-0.761,13.700},{-3.515,-3.013,13.940},{-2.775,-4.147,14.040},{-3.281,-5.104,14.163},{-1.416,-4.113,13.990},{-0.819,-5.022,14.061},{-8.189,-3.510,12.220},{-9.638,-3.338,12.000},{-7.536,-4.582,11.430},{-7.430,-2.130,11.970},{-8.600,1.610,10.390},{-5.841,0.470,10.610},{-5.602,1.448,10.192},{-7.039,-0.147,9.900},{-7.150,0.221,8.880},{-6.975,-1.234,9.851},{-8.195,0.300,10.780},{-9.046,-0.377,10.711},{-7.585,0.278,12.180},{-8.134,1.006,12.777},{-7.629,-1.044,12.910},{-6.843,-1.063,13.665},{-8.597,-1.152,13.399},{-6.186,0.650,11.970},{-1.223,-2.541,10.410},{-2.286,-1.912,10.460},{-2.298,-0.539,10.370},{-1.385,-0.076,10.262},{-3.408,0.280,10.410},{-3.321,1.492,10.320},{-4.615,-0.371,10.560},{-4.699,-1.738,10.660},{-5.674,-2.209,10.782},{-3.593,-2.506,10.610},{-3.627,-4.001,10.720},{-4.663,-4.341,10.734},{-3.112,-4.440,9.866},{-3.132,-4.309,11.640},{-8.688,1.973,8.840},{-9.759,2.965,8.620},{-8.790,0.723,8.050},{-7.263,2.644,8.590},{-6.012,6.358,7.010},{-4.449,3.814,7.230},{-3.681,4.464,6.812},{-5.781,4.018,6.520},{-5.655,4.382,5.500},{-6.369,3.101,6.471},{-6.453,5.060,7.400},{-7.540,5.013,7.331},{-5.973,4.683,8.800},{-5.990,5.595,9.397},{-6.786,3.639,9.530},{-6.160,3.163,10.285},{-7.632,4.121,10.019},{-4.622,4.162,8.590},{-2.483,-1.337,7.030},{-2.973,-0.203,7.080},{-2.176,0.914,6.990},{-1.165,0.753,6.882},{-2.593,2.230,7.030},{-1.809,3.159,6.940},{-3.952,2.412,7.180},{-4.823,1.356,7.280},{-5.889,1.548,7.402},{-4.379,0.084,7.230},{-5.286,-1.105,7.340},{-6.324,-0.771,7.354},{-5.128,-1.763,6.486},{-5.067,-1.646,8.260},{-5.869,6.703,5.460},{-6.153,8.135,5.240},{-6.686,5.751,4.670},{-4.322,6.408,5.210},{-1.126,8.677,3.630},{-1.358,5.700,3.850},{-0.354,5.775,3.432},{-2.315,6.649,3.140},{-1.999,6.869,2.120},{-3.329,6.253,3.091},{-2.247,7.887,4.020},{-3.154,8.487,3.951},{-2.080,7.300,5.420},{-1.557,8.047,6.017},{-3.351,6.932,6.150},{-3.125,6.180,6.905},{-3.752,7.820,6.639},{-1.293,6.084,5.210},{2.299,1.943,3.280},{2.894,2.784,3.273},{2.726,1.010,3.187},{-1.681,-0.315,3.510},{-1.143,0.792,3.530},{-1.735,2.068,3.660},{-3.077,2.430,3.800},{-3.046,3.735,3.880},{-3.943,4.343,3.999},{-1.779,4.275,3.800},{-0.934,3.196,3.660},{0.418,3.263,3.540},{0.974,2.061,3.410},{0.255,0.884,3.410},{0.783,0.005,3.315},{-0.808,8.873,2.080},{-0.196,10.198,1.860},{-2.029,8.583,1.290},{0.270,7.725,1.830},{4.189,7.682,0.250},{2.252,5.410,0.470},{3.108,4.880,0.052},{2.035,6.740,-0.240},{2.420,6.732,-1.260},{0.982,7.015,-0.289},{2.818,7.701,0.640},{2.437,8.720,0.571},{2.608,7.128,2.040},{3.470,7.426,2.637},{1.364,7.578,2.770},{1.104,6.836,3.525},{1.561,8.532,3.259},{2.530,5.682,1.830},{-2.039,1.948,0.270},{-1.112,2.765,0.320},{0.197,2.352,0.230},{0.356,1.340,0.122},{1.320,3.155,0.270},{2.445,2.697,0.180},{1.073,4.504,0.420},{-0.201,5.006,0.520},{-0.347,6.079,0.642},{-1.273,4.191,0.470},{-2.684,4.686,0.580},{-2.687,5.776,0.594},{-3.261,4.332,-0.274},{-3.131,4.310,1.500}};

MatrixForm[t1]

u={"P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","C","N","C","H","N","C","N","C","H","N","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","N","C","O","N","C","H","C","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","C","N","C","H","N","C","N","C","H","N","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","C","N","C","H","N","C","N","C","H","N","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","O","C","C","N","C","H","N","C","N","C","N","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","C","N","C","H","N","C","N","C","H","N","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","C","N","C","H","N","C","N","C","H","N","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","N","C","O","N","C","H","C","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","O","C","N","H","C","O","N","C","H","C","C","H","H","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","C","N","C","H","N","C","N","C","H","N","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","O","C","N","H","C","O","N","C","H","C","C","H","H","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","C","N","C","H","N","C","N","C","H","N","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","O","C","C","N","C","H","N","C","N","C","N","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","O","C","N","H","C","O","N","C","H","C","C","H","H","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","O","C","N","H","C","O","N","C","H","C","C","H","H","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","C","N","C","O","N","C","H","C","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","O","C","N","H","C","O","N","C","H","C","C","H","H","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","O","C","N","H","C","O","N","C","H","C","C","H","H","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","N","H","H","O","C","C","N","C","H","N","C","N","C","N","H","P","O","O","O","O","C","H","C","H","H","C","H","C","H","C","H","H","O","O","C","N","H","C","O","N","C","H","C","C","H","H","H"};

n=Length[t1]

un=Table[If[u[[i]]=="H","Z",u[[i]]],{i,1,n}]

num1=Table[i,{i,1,n}];

um=Table[{un[[i]],num1[[i]]},{i,1,n}];

us=Sort[um]

num=Table[us[[i,2]],{i,1,n}]

For[i=1,i<n+1,i++,If[us[[i,1]]=="Z",us[[i,1]]="H"]];

u1=Table[If[us[[i,1]]=="C",0.255,If[us[[i,1]]=="H",0.12,If[us[[i,1]]=="N",0.24,If[us[[i,1]]=="O",0.178,If[us[[i,1]]=="P",0.5]]]]],{i,1,n}]
u2=Table[If[us[[i,1]]=="C",0.255,If[us[[i,1]]=="H",0.6,If[us[[i,1]]=="N",0.22522,If[us[[i,1]]=="O",0.178,If[us[[i,1]]=="P",0.5407]]]]],{i,1,n}]

u=Table[us[[i,1]],{i,1,n}]

cn=Table[If[u[[i]]=="C",6.0,If[u[[i]]=="H",1.0,If[u[[i]]=="N",7.0,If[u[[i]]=="O",8.0,If[u[[i]]=="P",5.0]]]]],{i,1,n}]

{6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,6.`,2.,7.`,7.`,7.`,7.`,8.`,8.`,8.`,8.`,8.`,8.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`,1.`}

bc=0;
For[i=1, i<n+1,i++,If[u[[i]]≠"H", bc=bc+1]];
bc
bh=n-bc

DeleteFile[ "C:\\Kimb_Mathem\\tmp3"];
For[i=1,i<n+1,i++,
{u[[i]],Extract[t1,{i}]}>>> "C:\\Kimb_Mathem\\tmp3"]

t=Table[t1[[num[[i]]]],{i,1,n}];
Plot1o=Table[Graphics3D[{Opacity[0.8],If[u[[i]]=="N",Lighter[Blue],If[u[[i]]=="H",Lighter[White],If[u[[i]]=="O",Lighter[Red],If[u[[i]]=="Cl",Lighter[Green],If[u[[i]]=="P",Lighter[Magenta],Darker[Brown]]]]]],Sphere[{t[[i,1]],t[[i,2]],t[[i,3]]},u1[[i]]]}],{i,271,420}];
Plot2=Table[Graphics3D[Text[i,{t[[i,1]]-0.6,t[[i,2]],t[[i,3]]} ,{1,0}]],{i,271,420}];
Show[Plot1o,Plot2,{AspectRatio->Automatic,Boxed->True, Axes->Automatic,AxesLabel->{xl,yl,zl},SphericalRegion->True}]

(* Distance Matrix Bohr *)
dc=Table[Norm[t[[i]]-t[[j]]],{i,1,n},{j,1,n}];
MatrixForm[Prepend[dc,num1]]

x=Table[t[[i,1]],{i,1,n}];
y=Table[t[[i,2]],{i,1,n}];
z=Table[t[[i,3]],{i,1,n}];

(* connected pairs: distance criteria < 3.0 bohr *)
Dij=Table[If[(dc[[i,j]]<3.0||(u[[i]]=="P" && dc[[i,j]]<3.3)),dc[[i,j]],0],{i,1,n},{j,1,n}];
MatrixForm[Prepend[Dij,num1]]

(* count C atoms *)
nc=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="C",nc=nc+1,nc]]
cnt=Table[0,{nc}];
pri=cnt; sec=cnt;ter=cnt; qua=cnt;
nc

(* count O atoms *)
no=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="O",no=no+1,nc]]
no
(* count N atoms *)
nn=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="N",nn=nn+1,nc]]
nn
(* count P atoms *)
np=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="P",np=np+1,nc]]
np

(* count H atoms bonded to each C atom *)
For[j=1,j≤nc,j++,
For[i=bc+1,i≤n,i++,If[ Dij[[i,j]]≠0 ,cnt[[j]]=cnt[[j]]+1,0]]];
cnt

(* parse into prim, sec, tert, and quart C cores with different radii *)
n1=0; n2=0; n3=0; n4=0;
For[i=1,i≤nc,i++,
If[cnt[[i]]==3,n1=n1+1;pri[[n1]]=i,
If[cnt[[i]]==2,n2=n2+1;sec[[n2]]=i,
If[cnt[[i]]==1,n3=n3+1;ter[[n3]]=i,
If[cnt[[i]]=0,n4=n4+1;qua[[n4]]=i]]]]];
For [i=1, i≤nc,i++,If[pri[[i]]≠0,u2[[pri[[i]]]]=0.2572]]
For [i=1, i≤nc,i++,If[sec[[i]]≠0,u2[[sec[[i]]]]=0.2556]]
For [i=1, i≤nc,i++,If[ter[[i]]≠0,u2[[ter[[i]]]]=0.257]]
For [i=1, i≤nc,i++,If[qua[[i]]≠0,u2[[qua[[i]]]]=0.257]]

Eij=Table[If[i≥j,Dij[[i,j]],0],{i,1,nc},{j,1,nc}];
CD=Select[Flatten[Eij],(#≥2.4 &&#≤2.73) &]
Mean[CD] ± StandardDeviation[CD]

CC=Select[Flatten[Eij],(#≥2.73) &]
Mean[CC] ± StandardDeviation[CC]

Eij=Table[If[i≥j,Dij[[i,j]],0],{i,bc+1,n},{j,1,nc}];
CH=Select[Flatten[Eij],(#>0 &&#≤2.1) &]
Mean[CH] ± StandardDeviation[CH]

For[i=1 ,i<n+1,i++,
For[j=1, j<bc+1,j++,
If[(Dij[[i,j]]≠0  && u[[i]]== "H"),If[u[[j]]=="C",u2[[i]]=0.5693,If[u[[j]]=="F",u2[[i]]=0.5,If[u[[j]]=="N",u2[[i]]=0.5541,If[u[[j]]=="O",u2[[i]]=0.5241]]]]]]]
u2

(* subtract core radii from pair distances *)
Fij=Dij-u2;
(*MatrixForm[Fij];*)
Gij=Table[If[i≥j,If[Fij[[i,j]]>0,Fij[[i,j]],0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Gij];*)

Hij=Transpose[Gij];
(*MatrixForm[Hij];*)

Kij=Hij-u2;
(*MatrixForm[Kij];*)
Lij=Table[If[i≤j,If[Kij[[i,j]]>0,Kij[[i,j]],0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Lij];*)

(* determine bond cloud radii *)
Lij=Table[If[i≤j,If[Lij[[i,j]]>0,If[u[[j]]≠"H",0.5*Lij[[i,j]],Lij[[i,j]]],0]," "],{i,1,n},{j,1,n}];  
Qij=Lij;
(*MatrixForm[Lij];*)

(* determine kinetic energy *)
Nij=Table[If[i≤j,If[Lij[[i,j]]>0,2.25/Lij[[i,j]]^2,0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Nij]*)

Tval=Sum[If[i≤j,Nij[[i,j]],0],{i,1,n},{j,1,n}]

Tcor=Sum[If[u[[i]]≠"H",2.25/u2[[i]]^2,0],{i,1,n}]

Ttot=Tval+Tcor  (* plus π and lone pairs: + 8.7616 *)

TvalH=Sum[If[(i≤j && u[[j]]== "H"),Nij[[i,j]],0],{i,1,n},{j,1,n}]

TvalC=Tval-TvalH

Mij=Prepend[Qij,u];
Oij=Transpose[Mij];
Pij=Prepend[Oij,Prepend[u," "]];

con=Table[{i},{i,1,n}]; conj=Table[0,{nc}];
For[i=1,i≤n,i++,
For[j=1,j≤n,j++,If[ Dij[[i,j]]≠0 ,con[[i]]={con[[i]],j}]]];
conn=Table[Flatten[con[[i]]],{i,1,n}];
(* identify C atoms at double bonds *)
j=0;
For[i=1,i≤nc,i++,If[Length[conn[[i]]]==4,j=j+1;conj[[j]]=i]];
conj=Take[conj,{1,j}]
nj=Length[conj];
MatrixForm[conn, TableAlignments->Left]

Clear[c,v,p,Q,P,v1,v2,nk,nl]; (*c1=20; c2=21; c3=51;*) 
rpo=1.4;
rpi=1.4;
tpo=Table[{0,0,0},{nj+8}];
tpu=Table[{0,0,0},{nj+8}];
P=3; c={0,0,0,0,0}; 
(*--------------------------------------------------------------------------*)
(* construct the π clouds for C atoms *)

Do[(c=conn[[conj[[k]]]];
c=Take[c,P];
v=Table[Join[t[[c[[i]]]],{1}],{i,1,P}];
x2=t[[c[[2]],1]];
y2=t[[c[[2]],2]]; 
z2=t[[c[[2]],3]];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{-x2,-y2,-z2,1}};
p=v.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
T1=ArcTan[p[[1,3]]/p[[1,2]]];
U=p;
T=RotationMatrix[T1,{1,0,0}];
p=U.T;
T2=ArcTan[p[[1,2]]/p[[1,1]]];
U=p;
T=RotationMatrix[T2,{0,0,1}];
p=U.T;
T3=ArcTan[p[[3,3]]/p[[3,2]]];
U=p;
T=RotationMatrix[T3,{1,0,0}];
p=Chop[U.T];
mk=p;
ml=p;
mk[[1,3]]=rpo;
ml[[1,3]]=-rpi;
U=p;
U1=mk;
U2=ml;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
mk=U1.T;
ml=U2.T;
U=p;
U1=mk;
U2=ml;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
mk=U1.T;
ml=U2.T;
U=p;
U1=mk;
U2=ml;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
mk=U1.T;
ml=U2.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[mk[[i]],{1}],{i,1,P}];
v2=Table[Join[ml[[i]],{1}],{i,1,P}];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{x2,y2,z2,1}};
p=v.Q;
mk=v1.Q;
ml=v2.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
mk=Table[Take[mk[[i]],3],{i,1,P}];
ml=Table[Take[ml[[i]],3],{i,1,P}];
tpo[[k]]=mk[[1]];
tpu[[k]]=ml[[1]]),{k,nj}];

(*-------------------------------------------------------------------------*)
(* construct the π clouds for the 4 N atoms *)

k=1;
Clear[conn];
conn={{57,5,8},{58,13,16},{59,20,23},{60,27,30}};

Do[(c=conn[[k]];
v=Table[Join[t[[c[[i]]]],{1}],{i,1,P}];
x2=t[[c[[2]],1]];
y2=t[[c[[2]],2]]; 
z2=t[[c[[2]],3]];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{-x2,-y2,-z2,1}};
p=v.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
T1=ArcTan[p[[1,3]]/p[[1,2]]];
U=p;
T=RotationMatrix[T1,{1,0,0}];
p=U.T;
T2=ArcTan[p[[1,2]]/p[[1,1]]];
U=p;
T=RotationMatrix[T2,{0,0,1}];
p=U.T;
T3=ArcTan[p[[3,3]]/p[[3,2]]];
U=p;
T=RotationMatrix[T3,{1,0,0}];
p=Chop[U.T];
nk=p;
nl=p;
nk[[1,3]]=rpo;
nl[[1,3]]=-rpi;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[nk[[i]],{1}],{i,1,P}];
v2=Table[Join[nl[[i]],{1}],{i,1,P}];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{x2,y2,z2,1}};
p=v.Q;
nk=v1.Q;
nl=v2.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
nk=Table[Take[nk[[i]],3],{i,1,P}];
nl=Table[Take[nl[[i]],3],{i,1,P}];
tpo[[nj+k]]=nk[[1]];
tpu[[nj+k]]=nl[[1]]),{k,4}];

(*-------------------------------------------------------------------------*)
(* construct the π clouds for the 4 O atoms in C=O *)
k=1;
Clear[conn];
conn={{61,12,11},{63,18,15},{64,32,33},{66,34,35}};

Do[(c=conn[[k]];
v=Table[Join[t[[c[[i]]]],{1}],{i,1,P}];
x2=t[[c[[2]],1]];
y2=t[[c[[2]],2]]; 
z2=t[[c[[2]],3]];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{-x2,-y2,-z2,1}};
p=v.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
T1=ArcTan[p[[1,3]]/p[[1,2]]];
U=p;
T=RotationMatrix[T1,{1,0,0}];
p=U.T;
T2=ArcTan[p[[1,2]]/p[[1,1]]];
U=p;
T=RotationMatrix[T2,{0,0,1}];
p=U.T;
T3=ArcTan[p[[3,3]]/p[[3,2]]];
U=p;
T=RotationMatrix[T3,{1,0,0}];
p=Chop[U.T];
nk=p;
nl=p;
nk[[1,3]]=rpo;
nl[[1,3]]=-rpi;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[nk[[i]],{1}],{i,1,P}];
v2=Table[Join[nl[[i]],{1}],{i,1,P}];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{x2,y2,z2,1}};
p=v.Q;
nk=v1.Q;
nl=v2.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
nk=Table[Take[nk[[i]],3],{i,1,P}];
nl=Table[Take[nl[[i]],3],{i,1,P}];
tpo[[nj+4+k]]=nk[[1]];
tpu[[nj+4+k]]=nl[[1]]),{k,4}];

(* ---------------------------------------------------------------------------------------------
O atoms 62,66 LP *)

k=1;
Clear[conn];
conn={{36,62,12},{35,66,34}};
rlo=1.125;
tlo={{0.0,0.0,0.0},{0.0,0.0,0.0}};
tlu={{0.0,0.0,0.0},{0.0,0.0,0.0}};
Do[(c=conn[[k]];
a1=dc[[c[[2]],c[[1]]]]; a2=dc[[c[[2]],c[[3]]]]; d1=dc[[c[[1]],c[[3]]]];
α=ArcCos[(a1^2+a2^2-d1^2)/(2*a1*a2)];
L=(u2[[c[[2]]]]+rlo);
β=ArcSin[rlo/L];
If[k==1,γ=π/2-α/2,γ=3*π/2+α/2];
v=Table[Join[t[[c[[i]]]],{1}],{i,1,P}];
x2=t[[c[[2]],1]];
y2=t[[c[[2]],2]]; 
z2=t[[c[[2]],3]];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{-x2,-y2,-z2,1}};
p=v.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
T1=ArcTan[p[[1,3]]/p[[1,2]]];
U=p;
T=RotationMatrix[T1,{1,0,0}];
p=U.T;
T2=ArcTan[p[[1,2]]/p[[1,1]]];
U=p;
T=RotationMatrix[T2,{0,0,1}];
p=U.T;
T3=ArcTan[p[[3,3]]/p[[3,2]]];
U=p;
T=RotationMatrix[T3,{1,0,0}];
p=Chop[U.T];
nk=p;
nl=p;
nk[[2,2]]=L*Cos[γ]*Cos[β]; nk[[2,1]]=-L*Sin[γ]*Cos[β];  nk[[2,3]]=rlo; 
nl[[2,2]]=nk[[2,2]]; nl[[2,1]]=nk[[2,1]];  nl[[2,3]]=-rlo;
Plota = Table[Graphics3D[{Opacity[0.5],Sphere[p[[i]],u2[[c[[i]]]]]}],{i,1,P}];
Plotc = Table[Graphics3D[{Tube[{p[[i]],p[[j]]},0.02]}],{i,1,P-1},{j,i+1,P}];
Plotd=Graphics3D[{Opacity[0.5],Lighter[Red],Sphere[nk[[2]],rlo]}];
Plote=Graphics3D[{Opacity[0.5],Lighter[Red],Sphere[nl[[2]],rlo]}];
Plotb0=Table[Graphics3D[Text[c[[i]],p[[i]],{1,0}]],{i,1,P}];
Show[Plotd,Plote,Plota,Plotb0,(*Plotb1,Plotb2,*)Plotc,{Axes->Automatic,AxesLabel->{xl,yl,zl},AspectRatio->Automatic,Boxed->False,PlotRange->Automatic, SphericalRegion->True}]
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[nk[[i]],{1}],{i,1,P}];
v2=Table[Join[nl[[i]],{1}],{i,1,P}];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{x2,y2,z2,1}};
p=v.Q;
nk=v1.Q;
nl=v2.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
nk=Table[Take[nk[[i]],3],{i,1,P}];
nl=Table[Take[nl[[i]],3],{i,1,P}];
tlo[[k]]=nk[[2]];
tlu[[k]]=nl[[2]]),{k,2}];

(*----------------------------------------------------------------------------*)
(* next: O20 in C=O, two doubly occ LP in the same plane as the 3 atoms in c *)
c={4,20,18}
rlo=1.125;
a1=dc[[20,4]]; a2=dc[[20,18]]; d1=dc[[18,4]];
α=ArcCos[(a1^2+a2^2-d1^2)/(2*a1*a2)];
β=ArcSin[rlo/(u2[[20]]+rlo)];
γ=π/2-α/2;
L=(u2[[20]]+rlo)(**Cos[β]*);
(*c=Take[c,P];*)
v=Table[Join[t[[c[[i]]]],{1}],{i,1,P}];
x2=t[[c[[2]],1]];
y2=t[[c[[2]],2]]; 
z2=t[[c[[2]],3]];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{-x2,-y2,-z2,1}};
p=v.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
T1=ArcTan[p[[1,3]]/p[[1,2]]];
U=p;
T=RotationMatrix[T1,{1,0,0}];
p=U.T;
T2=ArcTan[p[[1,2]]/p[[1,1]]];
U=p;
T=RotationMatrix[T2,{0,0,1}];
p=U.T;
T3=ArcTan[p[[3,3]]/p[[3,2]]];
U=p;
T=RotationMatrix[T3,{1,0,0}];
p=Chop[U.T];
nk=p;
nl=p;
nk[[2,2]]=-L*Cos[γ+1.0472]; nk[[2,1]]=L*Sin[γ+1.0472];  nk[[2,3]]=0; 
nl[[2,2]]=-L*Cos[γ-1.0472]; nl[[2,1]]=L*Sin[γ-1.0472];  nl[[2,3]]=0;
(*Plota = Table[Graphics3D[{Opacity[0.5],Sphere[p[[i]],u2[[c[[i]]]]]}],{i,1,P}];
Plotc = Table[Graphics3D[{Tube[{p[[i]],p[[j]]},0.02]}],{i,1,P-1},{j,i+1,P}];
Plotd=Graphics3D[{Opacity[0.5],Lighter[Red],Sphere[nk[[2]],rlo]}];
Plote=Graphics3D[{Opacity[0.5],Lighter[Red],Sphere[nl[[2]],rlo]}];
Plotb0=Table[Graphics3D[Text[c[[i]],p[[i]],{1,0}]],{i,1,P}];
Show[Plotd,Plote,Plota,Plotb0,(*Plotb1,Plotb2,*)Plotc,{Axes->Automatic,AxesLabel->{xl,yl,zl},AspectRatio->Automatic,Boxed->False,PlotRange->All, SphericalRegion->True}]*)
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
nk=U1.T;
nl=U2.T;
U=p;
U1=nk;
U2=nl;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
nk=U1.T;
nl=U2.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[nk[[i]],{1}],{i,1,P}];
v2=Table[Join[nl[[i]],{1}],{i,1,P}];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{x2,y2,z2,1}};
p=v.Q;
nk=v1.Q;
nl=v2.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
nk=Table[Take[nk[[i]],3],{i,1,P}];
nl=Table[Take[nl[[i]],3],{i,1,P}];
tplo=nk[[2]];
tplu=nl[[2]];

(*-------------------------------------------------------------------------*)
(* next: N57, N58, N59, N60 one doubly occ LP as sp2 ligand *)
k=1;
Clear[conn];
tpln={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
conn={{5,57,8},{13,58,16},{20,59,23},{27,60,30}};
rln=1.246;
Do[(c=conn[[k]];
a1=dc[[c[[2]],c[[1]]]]; a2=dc[[c[[2]],c[[3]]]]; d1=dc[[c[[1]],c[[3]]]];
α=ArcCos[(a1^2+a2^2-d1^2)/(2*a1*a2)];
β=ArcSin[rln/(u2[[c[[2]]]]+rln)];
If[(k==1 || k==3),γ=π-α/4,γ=3*π/2+α/2];
L=(u2[[c[[2]]]]+rln);
v=Table[Join[t[[c[[i]]]],{1}],{i,1,P}];
x2=t[[c[[2]],1]];
y2=t[[c[[2]],2]]; 
z2=t[[c[[2]],3]];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{-x2,-y2,-z2,1}};
p=v.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
T1=ArcTan[p[[1,3]]/p[[1,2]]];
U=p;
T=RotationMatrix[T1,{1,0,0}];
p=U.T;
T2=ArcTan[p[[1,2]]/p[[1,1]]];
U=p;
T=RotationMatrix[T2,{0,0,1}];
p=U.T;
T3=ArcTan[p[[3,3]]/p[[3,2]]];
U=p;
T=RotationMatrix[T3,{1,0,0}];
p=Chop[U.T];
nk=p;
nk[[2,2]]=L*Cos[γ]; nk[[2,1]]=-L*Sin[γ];  nk[[2,3]]=0; 
(*Plota = Table[Graphics3D[{Opacity[0.5],Sphere[p[[i]],u2[[c[[i]]]]]}],{i,1,P}];
Plotc = Table[Graphics3D[{Tube[{p[[i]],p[[j]]},0.02]}],{i,1,P-1},{j,i+1,P}];
Plotd=Graphics3D[{Opacity[0.5],Lighter[Blue],Sphere[nk[[2]],rln]}];
Plotb0=Table[Graphics3D[Text[c[[i]],p[[i]],{1,0}]],{i,1,P}];
Show[Plotd,(*Plote,*)Plota,Plotb0,(*Plotb1,Plotb2,*)Plotc,{Axes->Automatic,AxesLabel->{xl,yl,zl},AspectRatio->Automatic,Boxed->False,PlotRange->All, SphericalRegion->True}]*)
U=p;
U1=nk;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
nk=U1.T;
U=p;
U1=nk;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
nk=U1.T;
U=p;
U1=nk;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
nk=U1.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[nk[[i]],{1}],{i,1,P}];
Q={{1,0,0,0},
{0,1,0,0},
{0,0,1,0},
{x2,y2,z2,1}};
p=v.Q;
nk=v1.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
nk=Table[Take[nk[[i]],3],{i,1,P}];
tpln[[k]]=nk[[2]]),{k,4}];

DeleteFile[ "C:\\Kimb_Mathem\\tmp4"];
For[i=2,i<n+2,i++,
For[j=2,j≤i,j++,
If[Pij[[i,j]]>0,{i-1,j-1,Extract[Pij,{i,j}]}>>> "C:\\Kimb_Mathem\\tmp4"]]]

a=ReadList["C:\\Kimb_Mathem\\tmp4"]
nb=Length[a]
cw=nb-bh

wxc=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(x[[a[[i,2]]]]-x[[a[[i,1]]]])+x[[a[[i,1]]]],{i,1,cw}];
wyc=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(y[[a[[i,2]]]]-y[[a[[i,1]]]])+y[[a[[i,1]]]],{i,1,cw}];   
wzc=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(z[[a[[i,2]]]]-z[[a[[i,1]]]])+z[[a[[i,1]]]],{i,1,cw}];

Plot7=Table[Graphics3D[{Darker[Brown,0.4],Opacity[0.5],If[wxc[[i]]≠0,Sphere[{wxc[[i]],wyc[[i]],wzc[[i]]},a[[i,3]]]]}],{i,1,cw}];

wxh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(x[[a[[i,1]]]]-x[[a[[i,2]]]])+x[[a[[i,2]]]],{i,cw+1,nb}];
wyh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(y[[a[[i,1]]]]-y[[a[[i,2]]]])+y[[a[[i,2]]]],{i,cw+1,nb}];
wzh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(z[[a[[i,1]]]]-z[[a[[i,2]]]])+z[[a[[i,2]]]],{i,cw+1,nb}];

Plot8=Table[Graphics3D[{Opacity[0.5],Sphere[{wxh[[i]],wyh[[i]],wzh[[i]]},a[[i+cw,3]]]}],{i,1,bh}];

Plot8a=Table[Graphics3D[{Opacity[0.3],Lighter[Cyan],Sphere[tpo[[i]],rpo]}],{i,1,nj+8}];
Plot8b=Table[Graphics3D[{Opacity[0.3],Lighter[Yellow],Sphere[tpu[[i]],rpi]}],{i,1,nj+8}];
Plot8c=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tlo[[i]],rlo]}],{i,1,2}];
Plot8d=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tlu[[i]],rlo]}],{i,1,2}];
Plot8e=Table[Graphics3D[{Opacity[0.4],Lighter[Blue],Sphere[tpln[[i]],rln]}],{i,1,4}];

Plot9= Table[Graphics3D[{Lighter[Orange,0.5],Tube[{{x[[a[[i,1]]]],y[[a[[i,1]]]],z[[a[[i,1]]]]},{x[[a[[i,2]]]],y[[a[[i,2]]]],z[[a[[i,2]]]]}},0.1]}],{i,1,nb}];
(*Plot9a= Graphics3D[Tube[{{x[[21]],y[[21]],z[[21]]},{nk[[2,1]],nk[[2,2]],nk[[2,3]]}},0.06]];
Plot9b= Graphics3D[Tube[{{x[[21]],y[[21]],z[[21]]},{nl[[2,1]],nl[[2,2]],nl[[2,3]]}},0.06]];*)

Show[Plot1o,(*Plot2,*)Plot7,Plot8,Plot8a,Plot8b,Plot8c,Plot8d,Plot8e,Plot9,{AspectRatio->Automatic,Boxed->False,(* Axes->Automatic,AxesLabel->{xl,yl,zl},*)SphericalRegion->True}]

Show[Plot1o,(*Plot2,*)Plot7,Plot8,Plot8a,Plot8b,Plot8c,Plot8d,Plot8e,Plot9,{AspectRatio->Automatic,Boxed->False,(* Axes->Automatic,AxesLabel->{xl,yl,zl},*)SphericalRegion->True}]

Show[Plot1o,Plot2,Plot9,(*Plot9a,Plot9b,*){Axes->Automatic,AxesLabel->{xl,yl,zl},AspectRatio->Automatic,Boxed->True, SphericalRegion->True}]

tc=Table[{x[[i]],y[[i]],z[[i]]},{i,1,bc}];   (* core clouds *)
twc=Table[{wxc[[i]],wyc[[i]],wzc[[i]]},{i,1,cw}];  (* Bond cöouds between cores *)
twh=Table[{wxh[[i]],wyh[[i]],wzh[[i]]},{i,1,bh}]; (* H-bonding clouds *) 
tw=Join[tc,twc,twh,tpo,tpu,tlo,tlu,tpln];
bt=Length[tw]
(*MatrixForm[tw]*)
(* coordinates of clouds *)

(* cloud-cloud distances *)
dw=Table[Norm[tw[[i]]-tw[[j]]],{i,1,bt},{j,1,bt}];
(*MatrixForm[dw]*)
Length[dw]

(* nucleus-nucleus distances *)
(*MatrixForm[dc] *)

tv1=Table[t[[i]],{i,1,bc}];
tv2=Table[{0,0,0},{i,1,cw}];
tv3= Table[t[[i]],{i,bc+1,n}];
tv4 =Table[{0,0,0},{i,1,nj+4}]; 
tv5 =Table[{0,0,0},{i,1,nj+4}]; 
tv6 =Table[{0,0,0},{i,1,16}]; 
tv  = Join[tv1,tv2,tv3,tv4,tv5,tv6];
MatrixForm[tv]
Length[tv]

ow=Join[Table[-2.0,{i,1,bc+nb}],Table[-0.5,{i,1,2*nj}],Table[-0.5,{i,1,6}],Table[-1.0,{i,1,4}],Table[-2.0,{i,1,16}]];
lw=Length[ow]
cnv=Join[Table[cn[[i]],{i,1,bc}],Table[0.0,{i,1,cw}],Table[cn[[i]],{i,bc+1,n}],Table[0.0,{i,1,2*(nj+5)}],Table[0.0,{i,1,16}]];
lc=Length[cnv]

echarg=Sum[ow[[i]],{i,1,lw}]
pcharg=Sum[cnv[[i]],{i,1,lc}]

(* nucleus-cloud distances *)
dcw=Table[sq=Norm[tw[[i]]-tv[[j]]];If[(sq> 0.1 && i≠j),1/sq,0],{i,1,bt},{j,1,bt}];
Length[dcw]
(*MatrixForm[dcw]*)

Vee=0.0; Vnn=0.0;

(* all clouds repulsion *)
For[i=1,i<bt,i++,
For[j=i+1,j<bt+1,j++,
Vee=Vee+ow[[i]]*ow[[j]]/dw[[i,j]]]]
Vee

(* half cloud with other common halfcloud: subtract WW *)
Vee=Vee-2/(rpo+rpi); (* 5*1/4 C + 3*1/4 N *)
Vee

sc=0.372;   (* core clouds *)
For[i=1 ,i<bc+1,i++,Vee=Vee+3*sc/u2[[i]]]
sw=0.378;  (* bonding clouds *)
For[i=1 ,i<nb+1,i++,Vee=Vee+3*sw/a[[i,3]]]
For[i=1 ,i<9,i++,Vee=Vee+3*sw/rlo]
For[i=1 ,i<4,i++,Vee=Vee+3*sw/rln]
Vee

(* all nuclei with all clouds, except i = j *)
Vne=cnv.(ow.dcw)

(* core nuclei in their cloud centers *)

For[i=1,i<bc+1,i++,Vne=Vne-3*cn[[i]]/u2[[i]]]
Vne

(* protons in their bonding cloud, distance u2[[a[[i,1]]]] to cloud center *)
For[i=cw+1,i<nb+1,i++,Vne=Vne-(3-(u2[[a[[i,1]]]]/a[[i,3]])^2)/a[[i,3]]]
Vne

(* nucleus-nucleus *)
For[i=1,i<n,i++,
For[j=i+1,j<n+1,j++,Vnn=Vnn+cn[[i]]*cn[[j]]/dc[[i,j]]]]
Vnn

Vee
Vnn
Vne
Tpi   = 2.296*(nj+3)*0.5*1.125*(1/rpo^2+1/rpi^2) (* pi kin e added *)
Tpi   = Tpi+ 2.296*24*1.0*1.125*(1/rpo^2+1/rpi^2) (* pi kin e added *)
Tlo   = 240*2*1.125*(1/rlo^2) (* lp kin e of 120 O added *)
Tln   = 73*2*1.125*(1/rln^2) (* lp kin e of 73 N added *)
Vtot=Vne+Vee+Vnn
Ttot=Ttot+Tpi+Tlo+Tln+20*337.06761
Vir   =-Vtot/(Ttot-20*337.06761)
Etot =Ttot+Vtot-40*337.06761