158
S. Rashev and D.C. Moule
F5566=-0.366164046715441, F5666=-0.358423904173009, F6666=0.0430217192246419;
double r_CO=1.2096, r_CH=1.1033; //equilibrium bond lengths
double Fact2=50341.17, x1, x2, x3, mt=0;
double z1=0, z2=0, z3=0, t1=0, t2=0, f=0;// input values for the calculation
// z1,z2 - C–H stretches; z3 - CO stretch; t1,t2 - in plane bending angles; f - o.p. angle
(stretches are in A, angles are in radians)
void main()
x1=z1/(z1+r_CH); x2=z2/(z2+r_CH); x3=z3/(z3+r_CO);
// harmonic
mt+=.5*F11*(x1+x2)*(x1+x2);
mt+=.5*F55*(x1-x2)*(x1-x2);
mt+=F22*x3*x3;
mt+=F44*f*f;
mt+=F12/sqrt(2.)*(x1+x2)*x3;
mt+=0.5*(F33+F66)*(t1*t1+t2*t2) + (F33-F66)*t1*t2;
mt+=F13/2.*(x1+x2)*(t1+t2);
mt+=F23/sqrt(2.)*x3*(t1+t2);
mt+=F56/2.*(x1-x2)*(t1-t2);
// cubic
mt+=F111/2./sqrt(2.)*(x1+x2)*(x1+x2)*(x1+x2);
mt+=F155/2./sqrt(2.)*(x1*x1-x2*x2)*(x1-x2);
mt+=F222*x3*x3*x3;
mt+=F112/2.*x3*(x1+x2)*(x1+x2);
mt+=F255/2.*x3*(x1-x2)*(x1-x2);
mt+=F122/sqrt(2.)*x3*x3*(x1+x2);
mt+=F144/sqrt(2.)*(x1+x2)*f*f;
mt+=F244*x3*f*f;
mt+=F113/2./sqrt(2.)*(x1+x2)*(x1+x2)*(t1+t2);
mt+=F355/2./sqrt(2.)*(x1-x2)*(x1-x2)*(t1+t2);
mt+=F133/2./sqrt(2.)*(x1+x2)*(t1+t2)*(t1+t2);
mt+=F156/2./sqrt(2.)*(x1*x1 - x2*x2)*(t1-t2);
mt+=F166/2./sqrt(2.)*(x1+x2)*(t1-t2)*(t1-t2);
mt+=F356/2./sqrt(2.)*(x1-x2)*(t1*t1 - t2*t2);
mt+=F223/sqrt(2.)*x3*x3*(t1+t2);
mt+=F233/2.*x3*(t1+t2)*(t1+t2);
mt+=F266/2.*x3*(t1-t2)*(t1-t2);
mt+=F123/2.*x3*(x1+x2)*(t1+t2);
mt+=F256/2.*x3*(x1-x2)*(t1-t2);
mt+=F344/sqrt(2.)*(t1+t2)*f*f;
mt+=F333/2./sqrt(2.)*(t1+t2)*(t1+t2)*(t1+t2);
mt+=F366/2./sqrt(2.)*(t1*t1 - t2*t2)*(t1-t2);
// quartic
mt+=F1111/4.*(x1+x2)*(x1+x2)*(x1+x2)*(x1+x2);
mt+=F5555/4.*(x1-x2)*(x1-x2)*(x1-x2)*(x1-x2);
mt+=F1155/4.*(x1*x1*x1*x1 - 2.*x1*x1*x2*x2 + x2*x2*x2*x2);
mt+=F1112/2./sqrt(2.)*x3*(x1+x2)*(x1+x2)*(x1+x2);