mengine/src/eopbend.c

#define EXTERN extern

#include "pcwin.h"
#include "pcmod.h"

#include "energies.h"
#include "angles.h"
#include "derivs.h"
#include "hess.h"
#include "utility.h"

void eopbend2b(int);

EXTERN struct t_minim_values {
        int iprint, ndc, nconst;
        float dielc;
        } minim_values;
        
void eopbend()
{
      int i,ia,ib,ic,id;
      double e,angle,dot,cosine;
      double dt,dt2,dt3,dt4;
      double xia,yia,zia;
      double xib,yib,zib;
      double xic,yic,zic;
      double xid,yid,zid;
      double xad,yad,zad;
      double xbd,ybd,zbd,rbd2;
      double xcd,ycd,zcd;
      double xpd,ypd,zpd,rpd2;
      double xt,yt,zt,rt2,delta;

      energies.eopb = 0.0;

      if (minim_values.iprint)
      {
          fprintf(pcmlogfile,"\nOut of Plane Bending Terms\n");
          fprintf(pcmlogfile,"    At1   At2   At3   At4      Angle     Oopconst     Eoop\n");
      }

      for (i=0; i < angles.nang; i++)
      {
          if (angles.angin[i] == TRUE)
          {
              ia = angles.i13[i][0];
              ib = angles.i13[i][1];
              ic = angles.i13[i][2];
              id = angles.i13[i][3];
              if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
              {
                 xia = atom[ia].x;
                 yia = atom[ia].y;
                 zia = atom[ia].z;
                 xib = atom[ib].x;
                 yib = atom[ib].y;
                 zib = atom[ib].z;
                 xic = atom[ic].x;
                 yic = atom[ic].y;
                 zic = atom[ic].z;
                 xid = atom[id].x;
                 yid = atom[id].y;
                 zid = atom[id].z;
                 xad = xia - xid;
                 yad = yia - yid;
                 zad = zia - zid;
                 xbd = xib - xid;
                 ybd = yib - yid;
                 zbd = zib - zid;
                 xcd = xic - xid;
                 ycd = yic - yid;
                 zcd = zic - zid;
                 xt = yad*zcd - zad*ycd;
                 yt = zad*xcd - xad*zcd;
                 zt = xad*ycd - yad*xcd;
                 rt2 = xt*xt + yt*yt + zt*zt;
                 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
                 xpd = xbd + xt*delta;
                 ypd = ybd + yt*delta;
                 zpd = zbd + zt*delta;
                 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
                 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
                 if (rbd2 != 0.0 && rpd2 != 0.0)
                 {
                    dot = xbd*xpd + ybd*ypd + zbd*zpd;
                    cosine = dot / sqrt(rbd2*rpd2);
                    if (cosine < -1.0)
                       cosine = -1.0;
                    if (cosine > 1.0)
                       cosine = 1.0;
                    angle = radian * acos(cosine);
                    dt = angle;
                    dt2 = dt * dt;
                    dt3 = dt2 * dt;
                    dt4 = dt2 * dt2;
                    e = units.angunit * angles.copb[i] * dt2
                        * (1.0+units.cang*dt+units.qang*dt2+units.pang*dt3+units.sang*dt4);
                    energies.eopb += e;
                    atom[ia].energy += e;
                    atom[ib].energy += e;
                    atom[ic].energy += e;
                    atom[id].energy += e;
                    if (minim_values.iprint)
                      fprintf(pcmlogfile,"Oop: %-4d- %-4d- %-4d- %-4d    %-8.3f    %-8.3f = %-8.4f\n",ia,ib,ic,id,
                        angle, angles.copb[i],e);
                 }
                 
              }
          }
      }
}

void eopbend1()
{
      int i,ia,ib,ic,id;
      double e,angle,dot,cosine;
      double dt,dt2,dt3,dt4;
      double deddt,termb,termp,term;
      double xia,yia,zia;
      double xib,yib,zib;
      double xic,yic,zic;
      double xid,yid,zid;
      double xad,yad,zad;
      double xbd,ybd,zbd,rbd2;
      double xcd,ycd,zcd;
      double xpd,ypd,zpd,rpd2;
      double xt,yt,zt,rt2,ptrt2;
      double xm,ym,zm,rm,delta,delta2;
      double dedxia,dedyia,dedzia;
      double dedxib,dedyib,dedzib;
      double dedxic,dedyic,dedzic;
      double dedxid,dedyid,dedzid;
      double dedxip,dedyip,dedzip;
      double dpdxia,dpdyia,dpdzia;
      double dpdxic,dpdyic,dpdzic;

      energies.eopb = 0.0;
      for (i=0; i <= natom; i++)
      {
          deriv.deopb[i][0] = 0.0;
          deriv.deopb[i][1] = 0.0;
          deriv.deopb[i][2] = 0.0;
      }

      for (i=0; i < angles.nang; i++)
      {
          if (angles.angin[i] == TRUE)
          {
              ia = angles.i13[i][0];
              ib = angles.i13[i][1];
              ic = angles.i13[i][2];
              id = angles.i13[i][3];
              if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
              {
                 xia = atom[ia].x;
                 yia = atom[ia].y;
                 zia = atom[ia].z;
                 xib = atom[ib].x;
                 yib = atom[ib].y;
                 zib = atom[ib].z;
                 xic = atom[ic].x;
                 yic = atom[ic].y;
                 zic = atom[ic].z;
                 xid = atom[id].x;
                 yid = atom[id].y;
                 zid = atom[id].z;
                 xad = xia - xid;
                 yad = yia - yid;
                 zad = zia - zid;
                 xbd = xib - xid;
                 ybd = yib - yid;
                 zbd = zib - zid;
                 xcd = xic - xid;
                 ycd = yic - yid;
                 zcd = zic - zid;
                 xt = yad*zcd - zad*ycd;
                 yt = zad*xcd - xad*zcd;
                 zt = xad*ycd - yad*xcd;
                 rt2 = xt*xt + yt*yt + zt*zt;
                 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
                 xpd = xbd + xt*delta;
                 ypd = ybd + yt*delta;
                 zpd = zbd + zt*delta;
                 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
                 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
                 xm = ypd*zbd - zpd*ybd;
                 ym = zpd*xbd - xpd*zbd;
                 zm = xpd*ybd - ypd*xbd;
                 rm = sqrt(xm*xm + ym*ym + zm*zm);
                 if (rm != 0.0)
                 {
                    dot = xbd*xpd + ybd*ypd + zbd*zpd;
                    cosine = dot / sqrt(rbd2*rpd2);
                    if (cosine < -1.0)
                       cosine = -1.0;
                    if (cosine > 1.0)
                       cosine = 1.0;
                    angle = radian * acos(cosine);
                    dt = angle;
                    dt2 = dt * dt;
                    dt3 = dt2 * dt;
                    dt4 = dt2 * dt2;
                    e = units.angunit * angles.copb[i] * dt2
                        * (1.0+units.cang*dt+units.qang*dt2+units.pang*dt3+units.sang*dt4);

                    deddt = units.angunit * angles.copb[i] * dt
                         * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
                             + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
                    deddt = deddt * radian;
                    termb = -deddt / (rbd2*rm);
                    termp = deddt / (rpd2*rm);
                    dedxib = termb * (ybd*zm-zbd*ym);
                    dedyib = termb * (zbd*xm-xbd*zm);
                    dedzib = termb * (xbd*ym-ybd*xm);
                    dedxip = termp * (ypd*zm-zpd*ym);
                    dedyip = termp * (zpd*xm-xpd*zm);
                    dedzip = termp * (xpd*ym-ypd*xm);

                    delta2 = 2.0 * delta;
                    ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
                    term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
                    dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
                    term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
                    dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
                    term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
                    dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
                    term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
                    dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
                    term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
                    dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
                    term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
                    dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;

                    dedxia = dpdxia;
                    dedyia = dpdyia;
                    dedzia = dpdzia;
                    dedxic = dpdxic;
                    dedyic = dpdyic;
                    dedzic = dpdzic;
                    dedxid = -dedxia - dedxib - dedxic;
                    dedyid = -dedyia - dedyib - dedyic;
                    dedzid = -dedzia - dedzib - dedzic;

                    energies.eopb += e;
                    deriv.deopb[ia][0] += dedxia;
                    deriv.deopb[ia][1] += dedyia;
                    deriv.deopb[ia][2] += dedzia;
                   
                    deriv.deopb[ib][0] += dedxib;
                    deriv.deopb[ib][1] += dedyib;
                    deriv.deopb[ib][2] += dedzib;
                    
                    deriv.deopb[ic][0] += dedxic;
                    deriv.deopb[ic][1] += dedyic;
                    deriv.deopb[ic][2] += dedzic;
                    
                    deriv.deopb[id][0] += dedxid;
                    deriv.deopb[id][1] += dedyid;
                    deriv.deopb[id][2] += dedzid;
                    virial.virx += xad*dedxia + xbd*dedxib + xcd*dedxic;
                    virial.viry += yad*dedyia + ybd*dedyib + ycd*dedyic;
                    virial.virz += zad*dedzia + zbd*dedzib + zcd*dedzic;
                 }
                 
              }
          }
      }
}

void eopbend2(int i)
{
    int j,k;
    int ia,ib,ic,id;
    double old,eps;
    double **d0; // d0[MAXATOM][3];

    d0 = dmatrix(0,natom,0,3);
    eps = 1.0e-7;

    for (k=0; k < angles.nang; k++)
    {
        if (angles.angin[k] == TRUE)
        {
            ia = angles.i13[k][0];
            ib = angles.i13[k][1];
            ic = angles.i13[k][2];
            id = angles.i13[k][3];
            if (ia == i || ib == i || ic == i || id == i)
            {
                eopbend2b(k);
                for (j = 0; j < 3; j++)
                {
                    d0[ia][j] = deriv.deopb[ia][j];
                    d0[ib][j] = deriv.deopb[ib][j];
                    d0[ic][j] = deriv.deopb[ic][j];
                    d0[id][j] = deriv.deopb[id][j];
                }
  // numerical x
                old = atom[i].x;
                atom[i].x += eps;
                eopbend2b(k);
                atom[i].x = old;
                for (j=0; j < 3; j++)
                {
                    hess.hessx[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
                    hess.hessx[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
                    hess.hessx[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
                    hess.hessx[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
                }
  // numerical y
                old = atom[i].y;
                atom[i].y += eps;
                eopbend2b(k);
                atom[i].y = old;
                for (j=0; j < 3; j++)
                {
                    hess.hessy[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
                    hess.hessy[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
                    hess.hessy[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
                    hess.hessy[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
                }
  // numerical z
                old = atom[i].z;
                atom[i].z += eps;
                eopbend2b(k);
                atom[i].z = old;
                for (j=0; j < 3; j++)
                {
                    hess.hessz[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
                    hess.hessz[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
                    hess.hessz[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
                    hess.hessz[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
                }
            }
        }
    }
    free_dmatrix( d0 ,0,natom,0,3);
}

void eopbend2b(int k)
{
    int j,ia,ib,ic,id;
    double angle,dot,cosine;
    double dt,dt2,dt3,dt4;
    double deddt,termb,termp,term;
    double xia,yia,zia;
    double xib,yib,zib;
    double xic,yic,zic;
    double xid,yid,zid;
    double xad,yad,zad;
    double xbd,ybd,zbd,rbd2;
    double xcd,ycd,zcd;
    double xpd,ypd,zpd,rpd2;
    double xt,yt,zt,rt2,ptrt2;
    double xm,ym,zm,rm,delta,delta2;
    double dedxia,dedyia,dedzia;
    double dedxib,dedyib,dedzib;
    double dedxic,dedyic,dedzic;
    double dedxid,dedyid,dedzid;
    double dedxip,dedyip,dedzip;
    double dpdxia,dpdyia,dpdzia;
    double dpdxic,dpdyic,dpdzic;
    
    ia = angles.i13[k][0];
    ib = angles.i13[k][1];
    ic = angles.i13[k][2];
    id = angles.i13[k][3];
                    
      xia = atom[ia].x;
      yia = atom[ia].y;
      zia = atom[ia].z;
      xib = atom[ib].x;
      yib = atom[ib].y;
      zib = atom[ib].z;
      xic = atom[ic].x;
      yic = atom[ic].y;
      zic = atom[ic].z;
      xid = atom[id].x;
      yid = atom[id].y;
      zid = atom[id].z;

      for (j=0; j < 3; j++)
      {
          deriv.deopb[ia][j] = 0.0;
          deriv.deopb[ib][j] = 0.0;
          deriv.deopb[ic][j] = 0.0;
          deriv.deopb[id][j] = 0.0;
      }
          
      xad = xia - xid;
      yad = yia - yid;
      zad = zia - zid;
      xbd = xib - xid;
      ybd = yib - yid;
      zbd = zib - zid;
      xcd = xic - xid;
      ycd = yic - yid;
      zcd = zic - zid;
      xt = yad*zcd - zad*ycd;
      yt = zad*xcd - xad*zcd;
      zt = xad*ycd - yad*xcd;
      rt2 = xt*xt + yt*yt + zt*zt;
      delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
      xpd = xbd + xt*delta;
      ypd = ybd + yt*delta;
      zpd = zbd + zt*delta;
      rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
      rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
      xm = ypd*zbd - zpd*ybd;
      ym = zpd*xbd - xpd*zbd;
      zm = xpd*ybd - ypd*xbd;
      rm = sqrt(xm*xm + ym*ym + zm*zm);
      if (rm != 0.0)
      {
         dot = xbd*xpd + ybd*ypd + zbd*zpd;
         cosine = dot / sqrt(rbd2*rpd2);
         if (cosine < -1.0)
            cosine = -1.0;
         if (cosine > 1.0)
            cosine = 1.0;
         angle = radian * acos(cosine);
         dt = angle;
         dt2 = dt * dt;
         dt3 = dt2 * dt;
         dt4 = dt2 * dt2;
         deddt = units.angunit * angles.copb[k] * dt
                   * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
                   + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
            
         deddt = deddt * radian;
         termb = -deddt / (rbd2*rm);
         termp = deddt / (rpd2*rm);
         dedxib = termb * (ybd*zm-zbd*ym);
         dedyib = termb * (zbd*xm-xbd*zm);
         dedzib = termb * (xbd*ym-ybd*xm);
         dedxip = termp * (ypd*zm-zpd*ym);
         dedyip = termp * (zpd*xm-xpd*zm);
         dedzip = termp * (xpd*ym-ypd*xm);

         delta2 = 2.0 * delta;
         ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
         term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
         dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
         term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
         dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
         term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
         dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
         term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
         dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
         term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
         dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
         term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
         dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;

         dedxia = dpdxia;
         dedyia = dpdyia;
         dedzia = dpdzia;
         dedxic = dpdxic;
         dedyic = dpdyic;
         dedzic = dpdzic;
         dedxid = -dedxia - dedxib - dedxic;
         dedyid = -dedyia - dedyib - dedyic;
         dedzid = -dedzia - dedzib - dedzic;

         deriv.deopb[ia][0] = dedxia;
         deriv.deopb[ia][1] = dedyia;
         deriv.deopb[ia][2] = dedzia;
         deriv.deopb[ib][0] = dedxib;
         deriv.deopb[ib][1] = dedyib;
         deriv.deopb[ib][2] = dedzib;
         deriv.deopb[ic][0] = dedxic;
         deriv.deopb[ic][1] = dedyic;
         deriv.deopb[ic][2] = dedzic;
         deriv.deopb[id][0] = dedxid;
         deriv.deopb[id][1] = dedyid;
         deriv.deopb[id][2] = dedzid;
      }
}
Revision:	93
Committed:	Sat Jan 17 23:24:55 2009 UTC (13 years, 4 months ago) by wdelano
File size:	16704 byte(s)
Log Message:	branch created
Line	File contents
1	#define EXTERN extern
2
3	#include "pcwin.h"
4	#include "pcmod.h"
5
6	#include "energies.h"
7	#include "angles.h"
8	#include "derivs.h"
9	#include "hess.h"
10	#include "utility.h"
11
12	void eopbend2b(int);
13
14	EXTERN struct t_minim_values {
15	int iprint, ndc, nconst;
16	float dielc;
17	} minim_values;
18
19	void eopbend()
20	{
21	int i,ia,ib,ic,id;
22	double e,angle,dot,cosine;
23	double dt,dt2,dt3,dt4;
24	double xia,yia,zia;
25	double xib,yib,zib;
26	double xic,yic,zic;
27	double xid,yid,zid;
28	double xad,yad,zad;
29	double xbd,ybd,zbd,rbd2;
30	double xcd,ycd,zcd;
31	double xpd,ypd,zpd,rpd2;
32	double xt,yt,zt,rt2,delta;
33
34	energies.eopb = 0.0;
35
36	if (minim_values.iprint)
37	{
38	fprintf(pcmlogfile,"\nOut of Plane Bending Terms\n");
39	fprintf(pcmlogfile," At1 At2 At3 At4 Angle Oopconst Eoop\n");
40	}
41
42	for (i=0; i < angles.nang; i++)
43	{
44	if (angles.angin[i] == TRUE)
45	{
46	ia = angles.i13[i][0];
47	ib = angles.i13[i][1];
48	ic = angles.i13[i][2];
49	id = angles.i13[i][3];
50	if (atom[ia].use \|\| atom[ib].use \|\| atom[ic].use \|\| atom[id].use)
51	{
52	xia = atom[ia].x;
53	yia = atom[ia].y;
54	zia = atom[ia].z;
55	xib = atom[ib].x;
56	yib = atom[ib].y;
57	zib = atom[ib].z;
58	xic = atom[ic].x;
59	yic = atom[ic].y;
60	zic = atom[ic].z;
61	xid = atom[id].x;
62	yid = atom[id].y;
63	zid = atom[id].z;
64	xad = xia - xid;
65	yad = yia - yid;
66	zad = zia - zid;
67	xbd = xib - xid;
68	ybd = yib - yid;
69	zbd = zib - zid;
70	xcd = xic - xid;
71	ycd = yic - yid;
72	zcd = zic - zid;
73	xt = yadzcd - zadycd;
74	yt = zadxcd - xadzcd;
75	zt = xadycd - yadxcd;
76	rt2 = xtxt + ytyt + zt*zt;
77	delta = -(xtxbd + ytybd + zt*zbd) / rt2;
78	xpd = xbd + xt*delta;
79	ypd = ybd + yt*delta;
80	zpd = zbd + zt*delta;
81	rbd2 = xbdxbd + ybdybd + zbd*zbd;
82	rpd2 = xpdxpd + ypdypd + zpd*zpd;
83	if (rbd2 != 0.0 && rpd2 != 0.0)
84	{
85	dot = xbdxpd + ybdypd + zbd*zpd;
86	cosine = dot / sqrt(rbd2*rpd2);
87	if (cosine < -1.0)
88	cosine = -1.0;
89	if (cosine > 1.0)
90	cosine = 1.0;
91	angle = radian * acos(cosine);
92	dt = angle;
93	dt2 = dt * dt;
94	dt3 = dt2 * dt;
95	dt4 = dt2 * dt2;
96	e = units.angunit * angles.copb[i] * dt2
97	* (1.0+units.cangdt+units.qangdt2+units.pangdt3+units.sangdt4);
98	energies.eopb += e;
99	atom[ia].energy += e;
100	atom[ib].energy += e;
101	atom[ic].energy += e;
102	atom[id].energy += e;
103	if (minim_values.iprint)
104	fprintf(pcmlogfile,"Oop: %-4d- %-4d- %-4d- %-4d %-8.3f %-8.3f = %-8.4f\n",ia,ib,ic,id,
105	angle, angles.copb[i],e);
106	}
107
108	}
109	}
110	}
111	}
112
113	void eopbend1()
114	{
115	int i,ia,ib,ic,id;
116	double e,angle,dot,cosine;
117	double dt,dt2,dt3,dt4;
118	double deddt,termb,termp,term;
119	double xia,yia,zia;
120	double xib,yib,zib;
121	double xic,yic,zic;
122	double xid,yid,zid;
123	double xad,yad,zad;
124	double xbd,ybd,zbd,rbd2;
125	double xcd,ycd,zcd;
126	double xpd,ypd,zpd,rpd2;
127	double xt,yt,zt,rt2,ptrt2;
128	double xm,ym,zm,rm,delta,delta2;
129	double dedxia,dedyia,dedzia;
130	double dedxib,dedyib,dedzib;
131	double dedxic,dedyic,dedzic;
132	double dedxid,dedyid,dedzid;
133	double dedxip,dedyip,dedzip;
134	double dpdxia,dpdyia,dpdzia;
135	double dpdxic,dpdyic,dpdzic;
136
137	energies.eopb = 0.0;
138	for (i=0; i <= natom; i++)
139	{
140	deriv.deopb[i][0] = 0.0;
141	deriv.deopb[i][1] = 0.0;
142	deriv.deopb[i][2] = 0.0;
143	}
144
145	for (i=0; i < angles.nang; i++)
146	{
147	if (angles.angin[i] == TRUE)
148	{
149	ia = angles.i13[i][0];
150	ib = angles.i13[i][1];
151	ic = angles.i13[i][2];
152	id = angles.i13[i][3];
153	if (atom[ia].use \|\| atom[ib].use \|\| atom[ic].use \|\| atom[id].use)
154	{
155	xia = atom[ia].x;
156	yia = atom[ia].y;
157	zia = atom[ia].z;
158	xib = atom[ib].x;
159	yib = atom[ib].y;
160	zib = atom[ib].z;
161	xic = atom[ic].x;
162	yic = atom[ic].y;
163	zic = atom[ic].z;
164	xid = atom[id].x;
165	yid = atom[id].y;
166	zid = atom[id].z;
167	xad = xia - xid;
168	yad = yia - yid;
169	zad = zia - zid;
170	xbd = xib - xid;
171	ybd = yib - yid;
172	zbd = zib - zid;
173	xcd = xic - xid;
174	ycd = yic - yid;
175	zcd = zic - zid;
176	xt = yadzcd - zadycd;
177	yt = zadxcd - xadzcd;
178	zt = xadycd - yadxcd;
179	rt2 = xtxt + ytyt + zt*zt;
180	delta = -(xtxbd + ytybd + zt*zbd) / rt2;
181	xpd = xbd + xt*delta;
182	ypd = ybd + yt*delta;
183	zpd = zbd + zt*delta;
184	rbd2 = xbdxbd + ybdybd + zbd*zbd;
185	rpd2 = xpdxpd + ypdypd + zpd*zpd;
186	xm = ypdzbd - zpdybd;
187	ym = zpdxbd - xpdzbd;
188	zm = xpdybd - ypdxbd;
189	rm = sqrt(xmxm + ymym + zm*zm);
190	if (rm != 0.0)
191	{
192	dot = xbdxpd + ybdypd + zbd*zpd;
193	cosine = dot / sqrt(rbd2*rpd2);
194	if (cosine < -1.0)
195	cosine = -1.0;
196	if (cosine > 1.0)
197	cosine = 1.0;
198	angle = radian * acos(cosine);
199	dt = angle;
200	dt2 = dt * dt;
201	dt3 = dt2 * dt;
202	dt4 = dt2 * dt2;
203	e = units.angunit * angles.copb[i] * dt2
204	* (1.0+units.cangdt+units.qangdt2+units.pangdt3+units.sangdt4);
205
206	deddt = units.angunit * angles.copb[i] * dt
207	* (2.0 + 3.0units.cangdt + 4.0units.qangdt2
208	+ 5.0units.pangdt3 + 6.0units.sangdt4);
209	deddt = deddt * radian;
210	termb = -deddt / (rbd2*rm);
211	termp = deddt / (rpd2*rm);
212	dedxib = termb * (ybdzm-zbdym);
213	dedyib = termb * (zbdxm-xbdzm);
214	dedzib = termb * (xbdym-ybdxm);
215	dedxip = termp * (ypdzm-zpdym);
216	dedyip = termp * (zpdxm-xpdzm);
217	dedzip = termp * (xpdym-ypdxm);
218
219	delta2 = 2.0 * delta;
220	ptrt2 = (dedxipxt + dedyipyt + dedzip*zt) / rt2;
221	term = (zcdybd-ycdzbd) + delta2(ytzcd-zt*ycd);
222	dpdxia = delta(ycddedzip-zcddedyip) + termptrt2;
223	term = (xcdzbd-zcdxbd) + delta2(ztxcd-xt*zcd);
224	dpdyia = delta(zcddedxip-xcddedzip) + termptrt2;
225	term = (ycdxbd-xcdybd) + delta2(xtycd-yt*xcd);
226	dpdzia = delta(xcddedyip-ycddedxip) + termptrt2;
227	term = (yadzbd-zadybd) + delta2(ztyad-yt*zad);
228	dpdxic = delta(zaddedyip-yaddedzip) + termptrt2;
229	term = (zadxbd-xadzbd) + delta2(xtzad-zt*xad);
230	dpdyic = delta(xaddedzip-zaddedxip) + termptrt2;
231	term = (xadybd-yadxbd) + delta2(ytxad-xt*yad);
232	dpdzic = delta(yaddedxip-xaddedyip) + termptrt2;
233
234	dedxia = dpdxia;
235	dedyia = dpdyia;
236	dedzia = dpdzia;
237	dedxic = dpdxic;
238	dedyic = dpdyic;
239	dedzic = dpdzic;
240	dedxid = -dedxia - dedxib - dedxic;
241	dedyid = -dedyia - dedyib - dedyic;
242	dedzid = -dedzia - dedzib - dedzic;
243
244	energies.eopb += e;
245	deriv.deopb[ia][0] += dedxia;
246	deriv.deopb[ia][1] += dedyia;
247	deriv.deopb[ia][2] += dedzia;
248
249	deriv.deopb[ib][0] += dedxib;
250	deriv.deopb[ib][1] += dedyib;
251	deriv.deopb[ib][2] += dedzib;
252
253	deriv.deopb[ic][0] += dedxic;
254	deriv.deopb[ic][1] += dedyic;
255	deriv.deopb[ic][2] += dedzic;
256
257	deriv.deopb[id][0] += dedxid;
258	deriv.deopb[id][1] += dedyid;
259	deriv.deopb[id][2] += dedzid;
260	virial.virx += xaddedxia + xbddedxib + xcd*dedxic;
261	virial.viry += yaddedyia + ybddedyib + ycd*dedyic;
262	virial.virz += zaddedzia + zbddedzib + zcd*dedzic;
263	}
264
265	}
266	}
267	}
268	}
269
270	void eopbend2(int i)
271	{
272	int j,k;
273	int ia,ib,ic,id;
274	double old,eps;
275	double **d0; // d0[MAXATOM][3];
276
277	d0 = dmatrix(0,natom,0,3);
278	eps = 1.0e-7;
279
280	for (k=0; k < angles.nang; k++)
281	{
282	if (angles.angin[k] == TRUE)
283	{
284	ia = angles.i13[k][0];
285	ib = angles.i13[k][1];
286	ic = angles.i13[k][2];
287	id = angles.i13[k][3];
288	if (ia == i \|\| ib == i \|\| ic == i \|\| id == i)
289	{
290	eopbend2b(k);
291	for (j = 0; j < 3; j++)
292	{
293	d0[ia][j] = deriv.deopb[ia][j];
294	d0[ib][j] = deriv.deopb[ib][j];
295	d0[ic][j] = deriv.deopb[ic][j];
296	d0[id][j] = deriv.deopb[id][j];
297	}
298	// numerical x
299	old = atom[i].x;
300	atom[i].x += eps;
301	eopbend2b(k);
302	atom[i].x = old;
303	for (j=0; j < 3; j++)
304	{
305	hess.hessx[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
306	hess.hessx[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
307	hess.hessx[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
308	hess.hessx[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
309	}
310	// numerical y
311	old = atom[i].y;
312	atom[i].y += eps;
313	eopbend2b(k);
314	atom[i].y = old;
315	for (j=0; j < 3; j++)
316	{
317	hess.hessy[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
318	hess.hessy[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
319	hess.hessy[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
320	hess.hessy[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
321	}
322	// numerical z
323	old = atom[i].z;
324	atom[i].z += eps;
325	eopbend2b(k);
326	atom[i].z = old;
327	for (j=0; j < 3; j++)
328	{
329	hess.hessz[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
330	hess.hessz[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
331	hess.hessz[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
332	hess.hessz[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
333	}
334	}
335	}
336	}
337	free_dmatrix( d0 ,0,natom,0,3);
338	}
339
340	void eopbend2b(int k)
341	{
342	int j,ia,ib,ic,id;
343	double angle,dot,cosine;
344	double dt,dt2,dt3,dt4;
345	double deddt,termb,termp,term;
346	double xia,yia,zia;
347	double xib,yib,zib;
348	double xic,yic,zic;
349	double xid,yid,zid;
350	double xad,yad,zad;
351	double xbd,ybd,zbd,rbd2;
352	double xcd,ycd,zcd;
353	double xpd,ypd,zpd,rpd2;
354	double xt,yt,zt,rt2,ptrt2;
355	double xm,ym,zm,rm,delta,delta2;
356	double dedxia,dedyia,dedzia;
357	double dedxib,dedyib,dedzib;
358	double dedxic,dedyic,dedzic;
359	double dedxid,dedyid,dedzid;
360	double dedxip,dedyip,dedzip;
361	double dpdxia,dpdyia,dpdzia;
362	double dpdxic,dpdyic,dpdzic;
363
364	ia = angles.i13[k][0];
365	ib = angles.i13[k][1];
366	ic = angles.i13[k][2];
367	id = angles.i13[k][3];
368
369	xia = atom[ia].x;
370	yia = atom[ia].y;
371	zia = atom[ia].z;
372	xib = atom[ib].x;
373	yib = atom[ib].y;
374	zib = atom[ib].z;
375	xic = atom[ic].x;
376	yic = atom[ic].y;
377	zic = atom[ic].z;
378	xid = atom[id].x;
379	yid = atom[id].y;
380	zid = atom[id].z;
381
382	for (j=0; j < 3; j++)
383	{
384	deriv.deopb[ia][j] = 0.0;
385	deriv.deopb[ib][j] = 0.0;
386	deriv.deopb[ic][j] = 0.0;
387	deriv.deopb[id][j] = 0.0;
388	}
389
390	xad = xia - xid;
391	yad = yia - yid;
392	zad = zia - zid;
393	xbd = xib - xid;
394	ybd = yib - yid;
395	zbd = zib - zid;
396	xcd = xic - xid;
397	ycd = yic - yid;
398	zcd = zic - zid;
399	xt = yadzcd - zadycd;
400	yt = zadxcd - xadzcd;
401	zt = xadycd - yadxcd;
402	rt2 = xtxt + ytyt + zt*zt;
403	delta = -(xtxbd + ytybd + zt*zbd) / rt2;
404	xpd = xbd + xt*delta;
405	ypd = ybd + yt*delta;
406	zpd = zbd + zt*delta;
407	rbd2 = xbdxbd + ybdybd + zbd*zbd;
408	rpd2 = xpdxpd + ypdypd + zpd*zpd;
409	xm = ypdzbd - zpdybd;
410	ym = zpdxbd - xpdzbd;
411	zm = xpdybd - ypdxbd;
412	rm = sqrt(xmxm + ymym + zm*zm);
413	if (rm != 0.0)
414	{
415	dot = xbdxpd + ybdypd + zbd*zpd;
416	cosine = dot / sqrt(rbd2*rpd2);
417	if (cosine < -1.0)
418	cosine = -1.0;
419	if (cosine > 1.0)
420	cosine = 1.0;
421	angle = radian * acos(cosine);
422	dt = angle;
423	dt2 = dt * dt;
424	dt3 = dt2 * dt;
425	dt4 = dt2 * dt2;
426	deddt = units.angunit * angles.copb[k] * dt
427	* (2.0 + 3.0units.cangdt + 4.0units.qangdt2
428	+ 5.0units.pangdt3 + 6.0units.sangdt4);
429
430	deddt = deddt * radian;
431	termb = -deddt / (rbd2*rm);
432	termp = deddt / (rpd2*rm);
433	dedxib = termb * (ybdzm-zbdym);
434	dedyib = termb * (zbdxm-xbdzm);
435	dedzib = termb * (xbdym-ybdxm);
436	dedxip = termp * (ypdzm-zpdym);
437	dedyip = termp * (zpdxm-xpdzm);
438	dedzip = termp * (xpdym-ypdxm);
439
440	delta2 = 2.0 * delta;
441	ptrt2 = (dedxipxt + dedyipyt + dedzip*zt) / rt2;
442	term = (zcdybd-ycdzbd) + delta2(ytzcd-zt*ycd);
443	dpdxia = delta(ycddedzip-zcddedyip) + termptrt2;
444	term = (xcdzbd-zcdxbd) + delta2(ztxcd-xt*zcd);
445	dpdyia = delta(zcddedxip-xcddedzip) + termptrt2;
446	term = (ycdxbd-xcdybd) + delta2(xtycd-yt*xcd);
447	dpdzia = delta(xcddedyip-ycddedxip) + termptrt2;
448	term = (yadzbd-zadybd) + delta2(ztyad-yt*zad);
449	dpdxic = delta(zaddedyip-yaddedzip) + termptrt2;
450	term = (zadxbd-xadzbd) + delta2(xtzad-zt*xad);
451	dpdyic = delta(xaddedzip-zaddedxip) + termptrt2;
452	term = (xadybd-yadxbd) + delta2(ytxad-xt*yad);
453	dpdzic = delta(yaddedxip-xaddedyip) + termptrt2;
454
455	dedxia = dpdxia;
456	dedyia = dpdyia;
457	dedzia = dpdzia;
458	dedxic = dpdxic;
459	dedyic = dpdyic;
460	dedzic = dpdzic;
461	dedxid = -dedxia - dedxib - dedxic;
462	dedyid = -dedyia - dedyib - dedyic;
463	dedzid = -dedzia - dedzib - dedzic;
464
465	deriv.deopb[ia][0] = dedxia;
466	deriv.deopb[ia][1] = dedyia;
467	deriv.deopb[ia][2] = dedzia;
468	deriv.deopb[ib][0] = dedxib;
469	deriv.deopb[ib][1] = dedyib;
470	deriv.deopb[ib][2] = dedzib;
471	deriv.deopb[ic][0] = dedxic;
472	deriv.deopb[ic][1] = dedyic;
473	deriv.deopb[ic][2] = dedzic;
474	deriv.deopb[id][0] = dedxid;
475	deriv.deopb[id][1] = dedyid;
476	deriv.deopb[id][2] = dedzid;
477	}
478	}