ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ebond.c
Revision: 103
Committed: Thu Feb 19 01:37:38 2009 UTC (10 years, 9 months ago) by gilbertke
File size: 5873 byte(s)
Log Message:
major rewrite - removing global data, adding electrostatics tag to read_sdf
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 (*)[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 void ebond2(int ia,int (*)[2],int (*)[MAXIAT],int (*)[MAXIAT],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[MAXBND][2],int *type,int *tclass,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, itype, ktype, classi,classk;
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 itype = type[it];
39 ktype = type[kt];
40 classi = tclass[it];
41 classk = tclass[kt];
42 bcorr = 0.0;
43 if ( use[it] || use[kt] )
44 {
45 xr = x[it] - x[kt];
46 yr = y[it] - y[kt];
47 zr = z[it] - z[kt];
48 rik2 = xr*xr + yr*yr + zr*zr;
49 rik = sqrt(rik2);
50 dt = rik - (bl[i] - bcorr);
51 dt2 = dt*dt;
52 e = units.bndunit *bk[i]*dt2*
53 (1.0 + units.cbnd*dt + units.qbnd*dt2);
54 *estr += e;
55 if (minim_values.iprint)
56 fprintf(pcmlogfile,"Bond: (%-3d) - (%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",it,kt,rik,bl[i]-bcorr,bk[i],e);
57 }
58 }
59 }
60 // =====================================
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 {
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 *estr = 0.0F;
70 for (i=0; i <= natom; i++)
71 {
72 deb[i][0] = 0.0;
73 deb[i][1] = 0.0;
74 deb[i][2] = 0.0;
75 }
76
77 for (i=0; i < nbnd; i++)
78 {
79 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 bcorr = 0.0;
86 if ( use[it] || use[kt] )
87 {
88 xr = x[it] - x[kt];
89 yr = y[it] - y[kt];
90 zr = z[it] - z[kt];
91 rik2 = xr*xr + yr*yr + zr*zr;
92 rik = sqrt(rik2);
93 dt = rik - (bl[i] - bcorr);
94 dt2 = dt*dt;
95 e = units.bndunit *bk[i]*dt2*
96 (1.0 + units.cbnd*dt + units.qbnd*dt2);
97 deddt = 2.0 * units.bndunit * bk[i] * dt
98 * (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 *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
114 }
115 }
116 }
117 // =============================================
118 void ebond2(int ia,int i12[MAXBND][2],int iat[MAXATOM][MAXIAT],int bo[MAXATOM][MAXIAT],double *x,double *y,double *z,double *bl,double *bk,float **hessx,float **hessy,float **hessz)
119 {
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 if (iat[ia][m]!= 0 && bo[ia][m] != 9)
128 {
129 ibond = find_bond(ia, iat[ia][m]);
130 if (i12[ibond][0] == ia)
131 k = i12[ibond][1];
132 else
133 k = i12[ibond][0];
134
135 xr = x[ia] - x[k];
136 yr = y[ia] - y[k];
137 zr = z[ia] - z[k];
138 rik2 = xr*xr + yr*yr + zr*zr;
139 rik = sqrt(rik2);
140 dt = rik - bl[ibond];
141 dt2 = dt * dt;
142
143 deddt = 2.0 * units.bndunit * bk[ibond] * dt
144 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
145 d2eddt2 = 2.0 * units.bndunit * bk[ibond]
146 * (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 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 }
180 }
181 }
182 }