ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ebond.c
Revision: 114
Committed: Wed Apr 29 19:07:31 2009 UTC (11 years, 2 months ago) by gilbertke
File size: 5428 byte(s)
Log Message:
added code for amber energy functions
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4
5 EXTERN struct t_minim_values {
6 int iprint, ndc, nconst;
7 float dielc;
8 } minim_values;
9
10 int find_bond(int,int);
11 // ebond(nbnd,i12,type,tclass,use,x,y,z,bl,bk,*estr)
12 // ebond1(nbnd,i12,type,tclass,use,x,y,z,bl,bk,*estr,**destr)
13 void ebond(int nbnd,int **,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr);
14 void ebond1(int natom,int nbnd,int **,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr,double **deb);
15 void ebond2(int ia,int **,int **,int **,double *x,double *y,double *z,double *bl,double *bk,float **hessx,float **hessy,float **hessz);
16
17 // ===========================
18 void ebond(int nbnd,int **i12,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr)
19 {
20 /* compute stretching energy */
21
22 int i, it, kt;
23 double xr, yr, zr, rik, rik2, bcorr;
24 double dt, dt2, e;
25
26 *estr = 0.0F;
27
28 if (minim_values.iprint)
29 {
30 fprintf(pcmlogfile,"\nBond Terms \n");
31 fprintf(pcmlogfile," At1 At2 R BLen Bconst Eb\n");
32 }
33
34 for (i=0; i < nbnd; i++)
35 {
36 it = i12[i][0];
37 kt = i12[i][1];
38 bcorr = 0.0;
39 if ( use[it] || use[kt] )
40 {
41 xr = x[it] - x[kt];
42 yr = y[it] - y[kt];
43 zr = z[it] - z[kt];
44 rik2 = xr*xr + yr*yr + zr*zr;
45 rik = sqrt(rik2);
46 dt = rik - (bl[i] - bcorr);
47 dt2 = dt*dt;
48 e = units.bndunit *bk[i]*dt2*
49 (1.0 + units.cbnd*dt + units.qbnd*dt2);
50 *estr += e;
51 if (minim_values.iprint)
52 fprintf(pcmlogfile,"Bond: (%-3d) - (%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",it,kt,rik,bl[i]-bcorr,bk[i],e);
53 }
54 }
55 }
56 // =====================================
57 void ebond1(int natom,int nbnd,int **i12,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr,double **deb)
58 {
59 /* compute stretching energy and first derivatives */
60 int i, it, kt;
61 double xr, yr, zr, rik, rik2, bcorr;
62 double dt, dt2, e, deddt;
63 double de,dedx,dedy,dedz;
64
65 *estr = 0.0F;
66 for (i=0; i <= natom; i++)
67 {
68 deb[i][0] = 0.0;
69 deb[i][1] = 0.0;
70 deb[i][2] = 0.0;
71 }
72
73 for (i=0; i < nbnd; i++)
74 {
75 it = i12[i][0];
76 kt = i12[i][1];
77 bcorr = 0.0;
78 if ( use[it] || use[kt] )
79 {
80 xr = x[it] - x[kt];
81 yr = y[it] - y[kt];
82 zr = z[it] - z[kt];
83 rik2 = xr*xr + yr*yr + zr*zr;
84 rik = sqrt(rik2);
85 dt = rik - (bl[i] - bcorr);
86 dt2 = dt*dt;
87 e = units.bndunit *bk[i]*dt2*
88 (1.0 + units.cbnd*dt + units.qbnd*dt2);
89 deddt = 2.0 * units.bndunit * bk[i] * dt
90 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
91 if (rik == 0.0)
92 de = 0.0;
93 else
94 de = deddt/rik;
95 dedx = de*xr;
96 dedy = de*yr;
97 dedz = de*zr;
98 *estr += e;
99 deb[it][0] += dedx;
100 deb[it][1] += dedy;
101 deb[it][2] += dedz;
102 deb[kt][0] -= dedx;
103 deb[kt][1] -= dedy;
104 deb[kt][2] -= dedz;
105
106 }
107 }
108 }
109 // =============================================
110 void ebond2(int ia,int **i12,int **iat,int **bo,double *x,double *y,double *z,double *bl,double *bk,float **hessx,float **hessy,float **hessz)
111 {
112 int j,k,m,ibond;
113 double dt,dt2;
114 double xr,yr,zr,rik,rik2,deddt,d2eddt2;
115 double de,term,termx,termy,termz,d2e[3][3];
116
117 for (m=0; m < MAXIAT; m++)
118 {
119 if (iat[ia][m]!= 0 && bo[ia][m] != 9)
120 {
121 ibond = find_bond(ia, iat[ia][m]);
122 if (i12[ibond][0] == ia)
123 k = i12[ibond][1];
124 else
125 k = i12[ibond][0];
126
127 xr = x[ia] - x[k];
128 yr = y[ia] - y[k];
129 zr = z[ia] - z[k];
130 rik2 = xr*xr + yr*yr + zr*zr;
131 rik = sqrt(rik2);
132 dt = rik - bl[ibond];
133 dt2 = dt * dt;
134
135 deddt = 2.0 * units.bndunit * bk[ibond] * dt
136 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
137 d2eddt2 = 2.0 * units.bndunit * bk[ibond]
138 * (1.0+3.0*units.cbnd*dt+6.0*units.qbnd*dt2);
139
140 if (rik2 == 0.0)
141 {
142 de = 0.0;
143 term = 0.0;
144 }else
145 {
146 de = deddt / rik;
147 term = (d2eddt2-de) / rik2;
148 }
149
150 termx = term * xr;
151 termy = term * yr;
152 termz = term * zr;
153 d2e[0][0] = termx*xr + de;
154 d2e[1][0] = termx*yr;
155 d2e[2][0] = termx*zr;
156 d2e[0][1] = d2e[1][0];
157 d2e[1][1] = termy*yr + de;
158 d2e[2][1] = termy*zr;
159 d2e[0][2] = d2e[2][0];
160 d2e[1][2] = d2e[2][1];
161 d2e[2][2] = termz*zr + de;
162
163 for (j=0; j < 3; j++)
164 {
165 hessx[ia][j] += d2e[j][0];
166 hessy[ia][j] += d2e[j][1];
167 hessz[ia][j] += d2e[j][2];
168 hessx[k][j] -= d2e[j][0];
169 hessy[k][j] -= d2e[j][1];
170 hessz[k][j] -= d2e[j][2];
171 }
172 }
173 }
174 }