ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ebond.c
Revision: 104
Committed: Fri Feb 20 14:09:46 2009 UTC (12 years, 7 months ago) by gilbertke
File size: 5825 byte(s)
Log Message:
full dynamic memory allocation of molecule
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     void ebond(int nbnd,int (*)[2],int *type,int *tclass,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr);
14     void ebond1(int natom,int nbnd,int (*)[2],int *type,int *tclass,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr,double **deb);
15 gilbertke 104 void ebond2(int ia,int (*)[2],int **,int **,double *x,double *y,double *z,double *bl,double *bk,float **hessx,float **hessy,float **hessz);
16 gilbertke 103
17     // ===========================
18     void ebond(int nbnd,int i12[MAXBND][2],int *type,int *tclass,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr)
19 wdelano 58 {
20     /* compute stretching energy */
21    
22     int i, it, kt, itype, ktype, classi,classk;
23     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     itype = type[it];
39     ktype = type[kt];
40     classi = tclass[it];
41     classk = tclass[kt];
42 wdelano 58 bcorr = 0.0;
43 gilbertke 103 if ( use[it] || use[kt] )
44 wdelano 58 {
45 gilbertke 103 xr = x[it] - x[kt];
46     yr = y[it] - y[kt];
47     zr = z[it] - z[kt];
48 wdelano 58 rik2 = xr*xr + yr*yr + zr*zr;
49     rik = sqrt(rik2);
50 gilbertke 103 dt = rik - (bl[i] - bcorr);
51 wdelano 58 dt2 = dt*dt;
52 gilbertke 103 e = units.bndunit *bk[i]*dt2*
53 wdelano 58 (1.0 + units.cbnd*dt + units.qbnd*dt2);
54 gilbertke 103 *estr += e;
55 wdelano 58 if (minim_values.iprint)
56 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);
57 wdelano 58 }
58     }
59     }
60 gilbertke 103 // =====================================
61     void ebond1(int natom,int nbnd,int i12[MAXBND][2],int *type,int *tclass,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr,double **deb)
62 wdelano 58 {
63     /* compute stretching energy and first derivatives */
64     int i, it, kt, itype, ktype, classi,classk;
65     double xr, yr, zr, rik, rik2, bcorr;
66     double dt, dt2, e, deddt;
67     double de,dedx,dedy,dedz;
68    
69 gilbertke 103 *estr = 0.0F;
70 wdelano 58 for (i=0; i <= natom; i++)
71     {
72 gilbertke 103 deb[i][0] = 0.0;
73     deb[i][1] = 0.0;
74     deb[i][2] = 0.0;
75 wdelano 58 }
76    
77 gilbertke 103 for (i=0; i < nbnd; i++)
78 wdelano 58 {
79 gilbertke 103 it = i12[i][0];
80     kt = i12[i][1];
81     itype = type[it];
82     ktype = type[kt];
83     classi = tclass[it];
84     classk = tclass[kt];
85 wdelano 58 bcorr = 0.0;
86 gilbertke 103 if ( use[it] || use[kt] )
87 wdelano 58 {
88 gilbertke 103 xr = x[it] - x[kt];
89     yr = y[it] - y[kt];
90     zr = z[it] - z[kt];
91 wdelano 58 rik2 = xr*xr + yr*yr + zr*zr;
92     rik = sqrt(rik2);
93 gilbertke 103 dt = rik - (bl[i] - bcorr);
94 wdelano 58 dt2 = dt*dt;
95 gilbertke 103 e = units.bndunit *bk[i]*dt2*
96 wdelano 58 (1.0 + units.cbnd*dt + units.qbnd*dt2);
97 gilbertke 103 deddt = 2.0 * units.bndunit * bk[i] * dt
98 wdelano 58 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
99     if (rik == 0.0)
100     de = 0.0;
101     else
102     de = deddt/rik;
103     dedx = de*xr;
104     dedy = de*yr;
105     dedz = de*zr;
106 gilbertke 103 *estr += e;
107     deb[it][0] += dedx;
108     deb[it][1] += dedy;
109     deb[it][2] += dedz;
110     deb[kt][0] -= dedx;
111     deb[kt][1] -= dedy;
112     deb[kt][2] -= dedz;
113 wdelano 58
114     }
115     }
116     }
117     // =============================================
118 gilbertke 104 void ebond2(int ia,int i12[MAXBND][2],int **iat,int **bo,double *x,double *y,double *z,double *bl,double *bk,float **hessx,float **hessy,float **hessz)
119 wdelano 58 {
120     int j,k,m,ibond;
121     double dt,dt2;
122     double xr,yr,zr,rik,rik2,deddt,d2eddt2;
123     double de,term,termx,termy,termz,d2e[3][3];
124    
125     for (m=0; m < MAXIAT; m++)
126     {
127 gilbertke 103 if (iat[ia][m]!= 0 && bo[ia][m] != 9)
128 wdelano 58 {
129 gilbertke 103 ibond = find_bond(ia, iat[ia][m]);
130     if (i12[ibond][0] == ia)
131     k = i12[ibond][1];
132 wdelano 58 else
133 gilbertke 103 k = i12[ibond][0];
134 wdelano 58
135 gilbertke 103 xr = x[ia] - x[k];
136     yr = y[ia] - y[k];
137     zr = z[ia] - z[k];
138 wdelano 58 rik2 = xr*xr + yr*yr + zr*zr;
139     rik = sqrt(rik2);
140 gilbertke 103 dt = rik - bl[ibond];
141 wdelano 58 dt2 = dt * dt;
142    
143 gilbertke 103 deddt = 2.0 * units.bndunit * bk[ibond] * dt
144 wdelano 58 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
145 gilbertke 103 d2eddt2 = 2.0 * units.bndunit * bk[ibond]
146 wdelano 58 * (1.0+3.0*units.cbnd*dt+6.0*units.qbnd*dt2);
147    
148     if (rik2 == 0.0)
149     {
150     de = 0.0;
151     term = 0.0;
152     }else
153     {
154     de = deddt / rik;
155     term = (d2eddt2-de) / rik2;
156     }
157    
158     termx = term * xr;
159     termy = term * yr;
160     termz = term * zr;
161     d2e[0][0] = termx*xr + de;
162     d2e[1][0] = termx*yr;
163     d2e[2][0] = termx*zr;
164     d2e[0][1] = d2e[1][0];
165     d2e[1][1] = termy*yr + de;
166     d2e[2][1] = termy*zr;
167     d2e[0][2] = d2e[2][0];
168     d2e[1][2] = d2e[2][1];
169     d2e[2][2] = termz*zr + de;
170    
171     for (j=0; j < 3; j++)
172     {
173 gilbertke 103 hessx[ia][j] += d2e[j][0];
174     hessy[ia][j] += d2e[j][1];
175     hessz[ia][j] += d2e[j][2];
176     hessx[k][j] -= d2e[j][0];
177     hessy[k][j] -= d2e[j][1];
178     hessz[k][j] -= d2e[j][2];
179 wdelano 58 }
180     }
181     }
182     }