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, 3 months ago) by gilbertke
File size: 5428 byte(s)
Log Message:
added code for amber energy functions
Line User Rev File contents
1 wdelano 58 #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 gilbertke 103 // 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 gilbertke 114 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 gilbertke 105 void ebond2(int ia,int **,int **,int **,double *x,double *y,double *z,double *bl,double *bk,float **hessx,float **hessy,float **hessz);
16 gilbertke 103
17     // ===========================
18 gilbertke 114 void ebond(int nbnd,int **i12,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr)
19 wdelano 58 {
20     /* compute stretching energy */
21    
22 gilbertke 114 int i, it, kt;
23 wdelano 58 double xr, yr, zr, rik, rik2, bcorr;
24     double dt, dt2, e;
25    
26 gilbertke 103 *estr = 0.0F;
27 wdelano 58
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 gilbertke 103 for (i=0; i < nbnd; i++)
35 wdelano 58 {
36 gilbertke 103 it = i12[i][0];
37     kt = i12[i][1];
38 wdelano 58 bcorr = 0.0;
39 gilbertke 103 if ( use[it] || use[kt] )
40 wdelano 58 {
41 gilbertke 103 xr = x[it] - x[kt];
42     yr = y[it] - y[kt];
43     zr = z[it] - z[kt];
44 wdelano 58 rik2 = xr*xr + yr*yr + zr*zr;
45     rik = sqrt(rik2);
46 gilbertke 103 dt = rik - (bl[i] - bcorr);
47 wdelano 58 dt2 = dt*dt;
48 gilbertke 103 e = units.bndunit *bk[i]*dt2*
49 wdelano 58 (1.0 + units.cbnd*dt + units.qbnd*dt2);
50 gilbertke 103 *estr += e;
51 wdelano 58 if (minim_values.iprint)
52 gilbertke 103 fprintf(pcmlogfile,"Bond: (%-3d) - (%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",it,kt,rik,bl[i]-bcorr,bk[i],e);
53 wdelano 58 }
54     }
55     }
56 gilbertke 103 // =====================================
57 gilbertke 114 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 wdelano 58 {
59     /* compute stretching energy and first derivatives */
60 gilbertke 114 int i, it, kt;
61 wdelano 58 double xr, yr, zr, rik, rik2, bcorr;
62     double dt, dt2, e, deddt;
63     double de,dedx,dedy,dedz;
64    
65 gilbertke 103 *estr = 0.0F;
66 wdelano 58 for (i=0; i <= natom; i++)
67     {
68 gilbertke 103 deb[i][0] = 0.0;
69     deb[i][1] = 0.0;
70     deb[i][2] = 0.0;
71 wdelano 58 }
72    
73 gilbertke 103 for (i=0; i < nbnd; i++)
74 wdelano 58 {
75 gilbertke 103 it = i12[i][0];
76     kt = i12[i][1];
77 wdelano 58 bcorr = 0.0;
78 gilbertke 103 if ( use[it] || use[kt] )
79 wdelano 58 {
80 gilbertke 103 xr = x[it] - x[kt];
81     yr = y[it] - y[kt];
82     zr = z[it] - z[kt];
83 wdelano 58 rik2 = xr*xr + yr*yr + zr*zr;
84     rik = sqrt(rik2);
85 gilbertke 103 dt = rik - (bl[i] - bcorr);
86 wdelano 58 dt2 = dt*dt;
87 gilbertke 103 e = units.bndunit *bk[i]*dt2*
88 wdelano 58 (1.0 + units.cbnd*dt + units.qbnd*dt2);
89 gilbertke 103 deddt = 2.0 * units.bndunit * bk[i] * dt
90 wdelano 58 * (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 gilbertke 103 *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 wdelano 58
106     }
107     }
108     }
109     // =============================================
110 gilbertke 105 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 wdelano 58 {
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 gilbertke 103 if (iat[ia][m]!= 0 && bo[ia][m] != 9)
120 wdelano 58 {
121 gilbertke 103 ibond = find_bond(ia, iat[ia][m]);
122     if (i12[ibond][0] == ia)
123     k = i12[ibond][1];
124 wdelano 58 else
125 gilbertke 103 k = i12[ibond][0];
126 wdelano 58
127 gilbertke 103 xr = x[ia] - x[k];
128     yr = y[ia] - y[k];
129     zr = z[ia] - z[k];
130 wdelano 58 rik2 = xr*xr + yr*yr + zr*zr;
131     rik = sqrt(rik2);
132 gilbertke 103 dt = rik - bl[ibond];
133 wdelano 58 dt2 = dt * dt;
134    
135 gilbertke 103 deddt = 2.0 * units.bndunit * bk[ibond] * dt
136 wdelano 58 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
137 gilbertke 103 d2eddt2 = 2.0 * units.bndunit * bk[ibond]
138 wdelano 58 * (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 gilbertke 103 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 wdelano 58 }
172     }
173     }
174     }