| 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 |
#include "field.h"
|
| 12 |
|
| 13 |
void eangle2a(int);
|
| 14 |
void eangle2b(int);
|
| 15 |
void eangle3(int, double (*)[3]);
|
| 16 |
//double acos(double);
|
| 17 |
|
| 18 |
|
| 19 |
EXTERN struct t_minim_values {
|
| 20 |
int iprint, ndc, nconst;
|
| 21 |
float dielc;
|
| 22 |
} minim_values;
|
| 23 |
|
| 24 |
void eangle()
|
| 25 |
{
|
| 26 |
int i, ia, ib, ic, id;
|
| 27 |
double e, angle, dot, cosine, factor;
|
| 28 |
double dt, dt2, dt3, dt4;
|
| 29 |
double xia, yia, zia;
|
| 30 |
double xib, yib, zib;
|
| 31 |
double xic, yic, zic;
|
| 32 |
double xid, yid, zid;
|
| 33 |
double xab, yab, zab, rab2;
|
| 34 |
double xcb, ycb, zcb, rcb2;
|
| 35 |
double xad, yad, zad;
|
| 36 |
double xbd, ybd, zbd;
|
| 37 |
double xcd,ycd,zcd;
|
| 38 |
double xip,yip,zip;
|
| 39 |
double xap,yap,zap,rap2;
|
| 40 |
double xcp,ycp,zcp,rcp2;
|
| 41 |
double xt,yt,zt,rt2,delta;
|
| 42 |
|
| 43 |
|
| 44 |
|
| 45 |
if (minim_values.iprint)
|
| 46 |
{
|
| 47 |
fprintf(pcmlogfile,"\nAngle Terms \n");
|
| 48 |
fprintf(pcmlogfile," At1 At2 At3 Angle Thet0 Tconst Ebend\n");
|
| 49 |
}
|
| 50 |
|
| 51 |
for (i=0; i < angles.nang; i++)
|
| 52 |
{
|
| 53 |
ia = angles.i13[i][0];
|
| 54 |
ib = angles.i13[i][1];
|
| 55 |
ic = angles.i13[i][2];
|
| 56 |
id = angles.i13[i][3];
|
| 57 |
if (atom[ia].use || atom[ib].use || atom[ic].use)
|
| 58 |
{
|
| 59 |
xia = atom[ia].x;
|
| 60 |
yia = atom[ia].y;
|
| 61 |
zia = atom[ia].z;
|
| 62 |
xib = atom[ib].x;
|
| 63 |
yib = atom[ib].y;
|
| 64 |
zib = atom[ib].z;
|
| 65 |
xic = atom[ic].x;
|
| 66 |
yic = atom[ic].y;
|
| 67 |
zic = atom[ic].z;
|
| 68 |
if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
|
| 69 |
{
|
| 70 |
xab = xia - xib;
|
| 71 |
yab = yia - yib;
|
| 72 |
zab = zia - zib;
|
| 73 |
rab2 = xab*xab + yab*yab + zab*zab;
|
| 74 |
xcb = xic - xib;
|
| 75 |
ycb = yic - yib;
|
| 76 |
zcb = zic - zib;
|
| 77 |
rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
|
| 78 |
if (rab2 != 0.0 && rcb2 != 0.00)
|
| 79 |
{
|
| 80 |
dot = xab*xcb + yab*ycb + zab*zcb;
|
| 81 |
cosine = dot / sqrt(rab2*rcb2);
|
| 82 |
if (cosine > 1.0)
|
| 83 |
cosine = 1.0;
|
| 84 |
if (cosine < -1.0)
|
| 85 |
cosine = -1.0;
|
| 86 |
angle = radian*acos(cosine);
|
| 87 |
|
| 88 |
if (angles.angtype[i] == HARMONIC)
|
| 89 |
{
|
| 90 |
dt = angle - angles.anat[i];
|
| 91 |
dt2 = dt*dt;
|
| 92 |
dt3 = dt2*dt;
|
| 93 |
dt4 = dt2*dt2;
|
| 94 |
e = units.angunit*angles.acon[i]*dt2
|
| 95 |
*(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
|
| 96 |
} else
|
| 97 |
{
|
| 98 |
factor = angle/radian;
|
| 99 |
e = units.angunit*angles.fcon[i]*
|
| 100 |
(angles.c0[i] + angles.c1[i]*cosine + angles.c2[i]*cos(2*factor)
|
| 101 |
+ angles.c3[i]*cos(3*factor) + angles.c4[i]*cos(4*factor)
|
| 102 |
+ angles.c5[i]*cos(5*factor) + angles.c6[i]*cos(6*factor));
|
| 103 |
}
|
| 104 |
energies.ebend += e;
|
| 105 |
atom[ia].energy += e;
|
| 106 |
atom[ib].energy += e;
|
| 107 |
atom[ic].energy += e;
|
| 108 |
if (minim_values.iprint)
|
| 109 |
{
|
| 110 |
if (angles.angtype[i] == HARMONIC)
|
| 111 |
fprintf(pcmlogfile,"Angle 1-3: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.4f = %-8.4f\n",
|
| 112 |
atom[ia].name,ia,atom[ib].name, ib,atom[ic].name, ic, angle, angles.anat[i],angles.acon[i], e);
|
| 113 |
else
|
| 114 |
fprintf(pcmlogfile,"Angle 1-3: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.4f = %-8.4f\n",
|
| 115 |
atom[ia].name,ia,atom[ib].name, ib,atom[ic].name, ic, angle, angles.fcon[i], e);
|
| 116 |
}
|
| 117 |
}
|
| 118 |
} else
|
| 119 |
{
|
| 120 |
xid = atom[id].x;
|
| 121 |
yid = atom[id].y;
|
| 122 |
zid = atom[id].z;
|
| 123 |
xad = xia - xid;
|
| 124 |
yad = yia - yid;
|
| 125 |
zad = zia - zid;
|
| 126 |
xbd = xib - xid;
|
| 127 |
ybd = yib - yid;
|
| 128 |
zbd = zib - zid;
|
| 129 |
xcd = xic - xid;
|
| 130 |
ycd = yic - yid;
|
| 131 |
zcd = zic - zid;
|
| 132 |
xt = yad*zcd - zad*ycd;
|
| 133 |
yt = zad*xcd - xad*zcd;
|
| 134 |
zt = xad*ycd - yad*xcd;
|
| 135 |
rt2 = xt*xt + yt*yt + zt*zt;
|
| 136 |
delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
|
| 137 |
xip = xib + xt*delta;
|
| 138 |
yip = yib + yt*delta;
|
| 139 |
zip = zib + zt*delta;
|
| 140 |
xap = xia - xip;
|
| 141 |
yap = yia - yip;
|
| 142 |
zap = zia - zip;
|
| 143 |
rap2 = xap*xap + yap*yap + zap*zap;
|
| 144 |
xcp = xic - xip;
|
| 145 |
ycp = yic - yip;
|
| 146 |
zcp = zic - zip;
|
| 147 |
rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
|
| 148 |
if (rap2 != 0.0 && rcp2 != 0.0)
|
| 149 |
{
|
| 150 |
dot = xap*xcp + yap*ycp + zap*zcp;
|
| 151 |
cosine = dot/ sqrt(rap2*rcp2);
|
| 152 |
if (cosine < -1.0)
|
| 153 |
cosine = -1.0;
|
| 154 |
if (cosine > 1.0)
|
| 155 |
cosine = 1.0;
|
| 156 |
angle = radian*acos(cosine);
|
| 157 |
dt = angle - angles.anat[i];
|
| 158 |
dt2 = dt*dt;
|
| 159 |
dt3 = dt2*dt;
|
| 160 |
dt4 = dt2*dt2;
|
| 161 |
e = units.angunit*angles.acon[i]*dt2
|
| 162 |
*(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
|
| 163 |
energies.ebend += e;
|
| 164 |
atom[ia].energy += e;
|
| 165 |
atom[ib].energy += e;
|
| 166 |
atom[ic].energy += e;
|
| 167 |
if (minim_values.iprint)
|
| 168 |
fprintf(pcmlogfile,"Angle InP: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.4f = %-8.4f\n",
|
| 169 |
atom[ia].name, ia, atom[ib].name, ib, atom[ic].name, ic, angle, angles.anat[i], angles.acon[i], e);
|
| 170 |
|
| 171 |
}
|
| 172 |
}
|
| 173 |
}
|
| 174 |
}
|
| 175 |
|
| 176 |
}
|
| 177 |
// ===============================================
|
| 178 |
void eangle1()
|
| 179 |
{
|
| 180 |
int i, ia, ib, ic, id;
|
| 181 |
double temp1, temp2;
|
| 182 |
double e, angle, dot, cosine, factor, sine;
|
| 183 |
double dt, dt2, dt3, dt4;
|
| 184 |
double deddt,terma,termc,term;
|
| 185 |
double xia, yia, zia;
|
| 186 |
double xib, yib, zib;
|
| 187 |
double xic, yic, zic;
|
| 188 |
double xid, yid, zid;
|
| 189 |
double xab, yab, zab, rab2;
|
| 190 |
double xcb, ycb, zcb, rcb2;
|
| 191 |
double xp,yp,zp,rp;
|
| 192 |
double xad, yad, zad;
|
| 193 |
double xbd, ybd, zbd;
|
| 194 |
double xcd,ycd,zcd;
|
| 195 |
double xip,yip,zip;
|
| 196 |
double xap,yap,zap,rap2;
|
| 197 |
double xcp,ycp,zcp,rcp2;
|
| 198 |
double xt,yt,zt,rt2,ptrt2;
|
| 199 |
double xm,ym,zm,rm,delta,delta2;
|
| 200 |
double dedxia,dedyia,dedzia;
|
| 201 |
double dedxib,dedyib,dedzib;
|
| 202 |
double dedxic,dedyic,dedzic;
|
| 203 |
double dedxid,dedyid,dedzid;
|
| 204 |
double dedxip,dedyip,dedzip;
|
| 205 |
double dpdxia,dpdyia,dpdzia;
|
| 206 |
double dpdxic,dpdyic,dpdzic;
|
| 207 |
double dtemp;
|
| 208 |
|
| 209 |
energies.ebend = 0.0F;
|
| 210 |
for (i=0; i <= natom; i++)
|
| 211 |
{
|
| 212 |
deriv.dea[i][0] = 0.0;
|
| 213 |
deriv.dea[i][1] = 0.0;
|
| 214 |
deriv.dea[i][2] = 0.0;
|
| 215 |
}
|
| 216 |
for (i=0; i < angles.nang; i++)
|
| 217 |
{
|
| 218 |
ia = angles.i13[i][0];
|
| 219 |
ib = angles.i13[i][1];
|
| 220 |
ic = angles.i13[i][2];
|
| 221 |
id = angles.i13[i][3];
|
| 222 |
if (atom[ia].use || atom[ib].use || atom[ic].use)
|
| 223 |
{
|
| 224 |
xia = atom[ia].x;
|
| 225 |
yia = atom[ia].y;
|
| 226 |
zia = atom[ia].z;
|
| 227 |
xib = atom[ib].x;
|
| 228 |
yib = atom[ib].y;
|
| 229 |
zib = atom[ib].z;
|
| 230 |
xic = atom[ic].x;
|
| 231 |
yic = atom[ic].y;
|
| 232 |
zic = atom[ic].z;
|
| 233 |
if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
|
| 234 |
{
|
| 235 |
xab = xia - xib;
|
| 236 |
yab = yia - yib;
|
| 237 |
zab = zia - zib;
|
| 238 |
rab2 = xab*xab + yab*yab + zab*zab;
|
| 239 |
xcb = xic - xib;
|
| 240 |
ycb = yic - yib;
|
| 241 |
zcb = zic - zib;
|
| 242 |
rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
|
| 243 |
xp = ycb*zab - zcb*yab;
|
| 244 |
yp = zcb*xab - xcb*zab;
|
| 245 |
zp = xcb*yab - ycb*xab;
|
| 246 |
rp = sqrt(xp*xp + yp*yp + zp*zp);
|
| 247 |
if (rp != 0.0 )
|
| 248 |
{
|
| 249 |
dot = xab*xcb + yab*ycb + zab*zcb;
|
| 250 |
cosine = dot / sqrt(rab2*rcb2);
|
| 251 |
if (cosine > 1.0)
|
| 252 |
cosine = 1.0;
|
| 253 |
if (cosine < -1.0)
|
| 254 |
cosine = -1.0;
|
| 255 |
angle = radian*acos(cosine);
|
| 256 |
if (angles.angtype[i] == HARMONIC)
|
| 257 |
{
|
| 258 |
dt = angle - angles.anat[i];
|
| 259 |
dt2 = dt*dt;
|
| 260 |
dt3 = dt2*dt;
|
| 261 |
dt4 = dt2*dt2;
|
| 262 |
|
| 263 |
e = units.angunit*angles.acon[i]*dt2
|
| 264 |
*(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
|
| 265 |
deddt = units.angunit*angles.acon[i]* dt * radian
|
| 266 |
* (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
|
| 267 |
+ 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
|
| 268 |
} else
|
| 269 |
{
|
| 270 |
factor = angle/radian;
|
| 271 |
sine = sin(factor);
|
| 272 |
e = units.angunit*angles.fcon[i]*
|
| 273 |
(angles.c0[i] + angles.c1[i]*cosine + angles.c2[i]*cos(2*factor)
|
| 274 |
+ angles.c3[i]*cos(3*factor) + angles.c4[i]*cos(4*factor)
|
| 275 |
+ angles.c5[i]*cos(5*factor) + angles.c6[i]*cos(6*factor));
|
| 276 |
deddt = -units.angunit*angles.fcon[i]*
|
| 277 |
(angles.c1[i]*sine + 2*angles.c2[i]*sin(2*factor)
|
| 278 |
+ 3*angles.c3[i]*sin(3*factor) + 4*angles.c4[i]*sin(4*factor)
|
| 279 |
+ 5*angles.c5[i]*sin(5*factor) + 6*angles.c6[i]*sin(6*factor));
|
| 280 |
}
|
| 281 |
|
| 282 |
/* terma = -deddt / (rab2*rp);
|
| 283 |
termc = deddt / (rcb2*rp);
|
| 284 |
dedxia = terma * (yab*zp-zab*yp);
|
| 285 |
dedyia = terma * (zab*xp-xab*zp);
|
| 286 |
dedzia = terma * (xab*yp-yab*xp);
|
| 287 |
dedxic = termc * (ycb*zp-zcb*yp);
|
| 288 |
dedyic = termc * (zcb*xp-xcb*zp);
|
| 289 |
dedzic = termc * (xcb*yp-ycb*xp);
|
| 290 |
dedxib = -dedxia - dedxic;
|
| 291 |
dedyib = -dedyia - dedyic;
|
| 292 |
dedzib = -dedzia - dedzic; */
|
| 293 |
|
| 294 |
dtemp = dot*dot/(rab2*rcb2);
|
| 295 |
if (fabs(dtemp) >= 1.0)
|
| 296 |
dtemp = 0.9999*dtemp/(fabs(dtemp));
|
| 297 |
terma = sqrt(1.0-dtemp);
|
| 298 |
termc = sqrt(rab2*rcb2);
|
| 299 |
dedxia = -deddt*( xcb/termc - (xab*cosine)/rab2)/terma;
|
| 300 |
dedyia = -deddt*( ycb/termc - (yab*cosine)/rab2)/terma;
|
| 301 |
dedzia = -deddt*( zcb/termc - (zab*cosine)/rab2)/terma;
|
| 302 |
|
| 303 |
dedxic = -deddt*( xab/termc - (xcb*cosine)/rcb2)/terma;
|
| 304 |
dedyic = -deddt*( yab/termc - (ycb*cosine)/rcb2)/terma;
|
| 305 |
dedzic = -deddt*( zab/termc - (zcb*cosine)/rcb2)/terma;
|
| 306 |
|
| 307 |
dedxib = -deddt*( (-xab-xcb)/termc + (dot*(xcb*rab2 + xab*rcb2))/(rab2*rcb2*termc))/terma;
|
| 308 |
dedyib = -deddt*( (-yab-ycb)/termc + (dot*(ycb*rab2 + yab*rcb2))/(rab2*rcb2*termc))/terma;
|
| 309 |
dedzib = -deddt*( (-zab-zcb)/termc + (dot*(zcb*rab2 + zab*rcb2))/(rab2*rcb2*termc))/terma;
|
| 310 |
|
| 311 |
|
| 312 |
deriv.dea[ia][0] += dedxia;
|
| 313 |
deriv.dea[ia][1] += dedyia;
|
| 314 |
deriv.dea[ia][2] += dedzia;
|
| 315 |
|
| 316 |
deriv.dea[ib][0] += dedxib;
|
| 317 |
deriv.dea[ib][1] += dedyib;
|
| 318 |
deriv.dea[ib][2] += dedzib;
|
| 319 |
|
| 320 |
deriv.dea[ic][0] += dedxic;
|
| 321 |
deriv.dea[ic][1] += dedyic;
|
| 322 |
deriv.dea[ic][2] += dedzic;
|
| 323 |
|
| 324 |
energies.ebend += e;
|
| 325 |
}
|
| 326 |
} else
|
| 327 |
{
|
| 328 |
xid = atom[id].x;
|
| 329 |
yid = atom[id].y;
|
| 330 |
zid = atom[id].z;
|
| 331 |
xad = xia - xid;
|
| 332 |
yad = yia - yid;
|
| 333 |
zad = zia - zid;
|
| 334 |
xbd = xib - xid;
|
| 335 |
ybd = yib - yid;
|
| 336 |
zbd = zib - zid;
|
| 337 |
xcd = xic - xid;
|
| 338 |
ycd = yic - yid;
|
| 339 |
zcd = zic - zid;
|
| 340 |
xt = yad*zcd - zad*ycd;
|
| 341 |
yt = zad*xcd - xad*zcd;
|
| 342 |
zt = xad*ycd - yad*xcd;
|
| 343 |
rt2 = xt*xt + yt*yt + zt*zt;
|
| 344 |
delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
|
| 345 |
xip = xib + xt*delta;
|
| 346 |
yip = yib + yt*delta;
|
| 347 |
zip = zib + zt*delta;
|
| 348 |
xap = xia - xip;
|
| 349 |
yap = yia - yip;
|
| 350 |
zap = zia - zip;
|
| 351 |
rap2 = xap*xap + yap*yap + zap*zap;
|
| 352 |
xcp = xic - xip;
|
| 353 |
ycp = yic - yip;
|
| 354 |
zcp = zic - zip;
|
| 355 |
rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
|
| 356 |
xm = ycp*zap - zcp*yap;
|
| 357 |
ym = zcp*xap - xcp*zap;
|
| 358 |
zm = xcp*yap - ycp*xap;
|
| 359 |
rm = sqrt(xm*xm + ym*ym + zm*zm);
|
| 360 |
if (rm != 0.0)
|
| 361 |
{
|
| 362 |
dot = xap*xcp + yap*ycp + zap*zcp;
|
| 363 |
cosine = dot/ sqrt(rap2*rcp2);
|
| 364 |
if (cosine < -1.0)
|
| 365 |
cosine = 1.0;
|
| 366 |
if (cosine > 1.0)
|
| 367 |
cosine = 1.0;
|
| 368 |
temp1 = acos(cosine);
|
| 369 |
temp2 = radian;
|
| 370 |
angle = temp1*temp2;
|
| 371 |
dt = angle - angles.anat[i];
|
| 372 |
dt2 = dt*dt;
|
| 373 |
dt3 = dt2*dt;
|
| 374 |
dt4 = dt2*dt2;
|
| 375 |
e = units.angunit*angles.acon[i]*dt2
|
| 376 |
*(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
|
| 377 |
|
| 378 |
deddt = units.angunit *angles.acon[i]* dt
|
| 379 |
* (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
|
| 380 |
+ 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
|
| 381 |
|
| 382 |
deddt = deddt * radian;
|
| 383 |
terma = -deddt / (rap2*rm);
|
| 384 |
termc = deddt / (rcp2*rm);
|
| 385 |
dedxia = terma * (yap*zm-zap*ym);
|
| 386 |
dedyia = terma * (zap*xm-xap*zm);
|
| 387 |
dedzia = terma * (xap*ym-yap*xm);
|
| 388 |
dedxic = termc * (ycp*zm-zcp*ym);
|
| 389 |
dedyic = termc * (zcp*xm-xcp*zm);
|
| 390 |
dedzic = termc * (xcp*ym-ycp*xm);
|
| 391 |
dedxip = -dedxia - dedxic;
|
| 392 |
dedyip = -dedyia - dedyic;
|
| 393 |
dedzip = -dedzia - dedzic;
|
| 394 |
|
| 395 |
delta2 = 2.0 * delta;
|
| 396 |
ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
|
| 397 |
term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
|
| 398 |
dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
|
| 399 |
term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
|
| 400 |
dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
|
| 401 |
term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
|
| 402 |
dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
|
| 403 |
term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
|
| 404 |
dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
|
| 405 |
term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
|
| 406 |
dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
|
| 407 |
term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
|
| 408 |
dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
|
| 409 |
|
| 410 |
dedxia = dedxia + dpdxia;
|
| 411 |
dedyia = dedyia + dpdyia;
|
| 412 |
dedzia = dedzia + dpdzia;
|
| 413 |
dedxib = dedxip;
|
| 414 |
dedyib = dedyip;
|
| 415 |
dedzib = dedzip;
|
| 416 |
dedxic = dedxic + dpdxic;
|
| 417 |
dedyic = dedyic + dpdyic;
|
| 418 |
dedzic = dedzic + dpdzic;
|
| 419 |
dedxid = -dedxia - dedxib - dedxic;
|
| 420 |
dedyid = -dedyia - dedyib - dedyic;
|
| 421 |
dedzid = -dedzia - dedzib - dedzic;
|
| 422 |
|
| 423 |
energies.ebend += e;
|
| 424 |
deriv.dea[ia][0] += dedxia;
|
| 425 |
deriv.dea[ia][1] += dedyia;
|
| 426 |
deriv.dea[ia][2] += dedzia;
|
| 427 |
|
| 428 |
deriv.dea[ib][0] += dedxib;
|
| 429 |
deriv.dea[ib][1] += dedyib;
|
| 430 |
deriv.dea[ib][2] += dedzib;
|
| 431 |
|
| 432 |
deriv.dea[ic][0] += dedxic;
|
| 433 |
deriv.dea[ic][1] += dedyic;
|
| 434 |
deriv.dea[ic][2] += dedzic;
|
| 435 |
|
| 436 |
deriv.dea[id][0] += dedxid;
|
| 437 |
deriv.dea[id][1] += dedyid;
|
| 438 |
deriv.dea[id][2] += dedzid;
|
| 439 |
|
| 440 |
}
|
| 441 |
}
|
| 442 |
}
|
| 443 |
}
|
| 444 |
}
|
| 445 |
|
| 446 |
|
| 447 |
void eangle2(int i)
|
| 448 |
{
|
| 449 |
int j,k,ia,ib,ic,id;
|
| 450 |
double old,eps;
|
| 451 |
double **d0; // d0[MAXATOM][3];
|
| 452 |
|
| 453 |
d0 = dmatrix(0,natom,0,3);
|
| 454 |
|
| 455 |
eangle2a(i);
|
| 456 |
|
| 457 |
eps = 1.0e-7;
|
| 458 |
for (k=0; k < angles.nang; k++)
|
| 459 |
{
|
| 460 |
if ( (angles.angin[k] == TRUE) && (field.type != MMFF94) )
|
| 461 |
{
|
| 462 |
ia = angles.i13[k][0];
|
| 463 |
ib = angles.i13[k][1];
|
| 464 |
ic = angles.i13[k][2];
|
| 465 |
id = angles.i13[k][3];
|
| 466 |
if (i == ia || i == ib || i == ic || i == id)
|
| 467 |
{
|
| 468 |
eangle2b(k);
|
| 469 |
for (j=0; j < 3; j++)
|
| 470 |
{
|
| 471 |
d0[ia][j] = deriv.dea[ia][j];
|
| 472 |
d0[ib][j] = deriv.dea[ib][j];
|
| 473 |
d0[ic][j] = deriv.dea[ic][j];
|
| 474 |
d0[id][j] = deriv.dea[id][j];
|
| 475 |
}
|
| 476 |
// numerical x derivative
|
| 477 |
old = atom[i].x;
|
| 478 |
atom[i].x += eps;
|
| 479 |
eangle2b(k);
|
| 480 |
atom[i].x = old;
|
| 481 |
for (j=0; j < 3; j++)
|
| 482 |
{
|
| 483 |
hess.hessx[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
|
| 484 |
hess.hessx[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
|
| 485 |
hess.hessx[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
|
| 486 |
hess.hessx[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
|
| 487 |
}
|
| 488 |
// numerical y derivative
|
| 489 |
old = atom[i].y;
|
| 490 |
atom[i].y += eps;
|
| 491 |
eangle2b(k);
|
| 492 |
atom[i].y = old;
|
| 493 |
for (j=0; j < 3; j++)
|
| 494 |
{
|
| 495 |
hess.hessy[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
|
| 496 |
hess.hessy[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
|
| 497 |
hess.hessy[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
|
| 498 |
hess.hessy[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
|
| 499 |
}
|
| 500 |
// numerical z derivative
|
| 501 |
old = atom[i].z;
|
| 502 |
atom[i].z += eps;
|
| 503 |
eangle2b(k);
|
| 504 |
atom[i].z = old;
|
| 505 |
for (j=0; j < 3; j++)
|
| 506 |
{
|
| 507 |
hess.hessz[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
|
| 508 |
hess.hessz[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
|
| 509 |
hess.hessz[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
|
| 510 |
hess.hessz[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
|
| 511 |
}
|
| 512 |
}
|
| 513 |
}
|
| 514 |
}
|
| 515 |
free_dmatrix(d0 ,0,natom,0,3);
|
| 516 |
}
|
| 517 |
// =================================================
|
| 518 |
void eangle2a(int iatom)
|
| 519 |
{
|
| 520 |
int i,ia,ib,ic;
|
| 521 |
double angle,dot,cosine,factor,sine;
|
| 522 |
double dt,dt2,dt3,dt4;
|
| 523 |
double xia,yia,zia;
|
| 524 |
double xib,yib,zib;
|
| 525 |
double xic,yic,zic;
|
| 526 |
double xab,yab,zab,rab;
|
| 527 |
double xcb,ycb,zcb,rcb;
|
| 528 |
double xp,yp,zp,rp,rp2;
|
| 529 |
double xrab,yrab,zrab,rab2;
|
| 530 |
double xrcb,yrcb,zrcb,rcb2;
|
| 531 |
double xabp,yabp,zabp;
|
| 532 |
double xcbp,ycbp,zcbp;
|
| 533 |
double deddt,d2eddt2,terma,termc;
|
| 534 |
double ddtdxia,ddtdyia,ddtdzia;
|
| 535 |
double ddtdxib,ddtdyib,ddtdzib;
|
| 536 |
double ddtdxic,ddtdyic,ddtdzic;
|
| 537 |
double dxiaxia,dxiayia,dxiazia;
|
| 538 |
double dxibxib,dxibyib,dxibzib;
|
| 539 |
double dxicxic,dxicyic,dxiczic;
|
| 540 |
double dyiayia,dyiazia,dziazia;
|
| 541 |
double dyibyib,dyibzib,dzibzib;
|
| 542 |
double dyicyic,dyiczic,dziczic;
|
| 543 |
double dxibxia,dxibyia,dxibzia;
|
| 544 |
double dyibxia,dyibyia,dyibzia;
|
| 545 |
double dzibxia,dzibyia,dzibzia;
|
| 546 |
double dxibxic,dxibyic,dxibzic;
|
| 547 |
double dyibxic,dyibyic,dyibzic;
|
| 548 |
double dzibxic,dzibyic,dzibzic;
|
| 549 |
double dxiaxic,dxiayic,dxiazic;
|
| 550 |
double dyiaxic,dyiayic,dyiazic;
|
| 551 |
double dziaxic,dziayic,dziazic;
|
| 552 |
|
| 553 |
for (i=0; i < angles.nang; i++)
|
| 554 |
{
|
| 555 |
ia = angles.i13[i][0];
|
| 556 |
ib = angles.i13[i][1];
|
| 557 |
ic = angles.i13[i][2];
|
| 558 |
if (iatom == ia || iatom == ib || iatom == ic)
|
| 559 |
{
|
| 560 |
xia = atom[ia].x;
|
| 561 |
yia = atom[ia].y;
|
| 562 |
zia = atom[ia].z;
|
| 563 |
xib = atom[ib].x;
|
| 564 |
yib = atom[ib].y;
|
| 565 |
zib = atom[ib].z;
|
| 566 |
xic = atom[ic].x;
|
| 567 |
yic = atom[ic].y;
|
| 568 |
zic = atom[ic].z;
|
| 569 |
|
| 570 |
if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
|
| 571 |
{
|
| 572 |
xab = xia - xib;
|
| 573 |
yab = yia - yib;
|
| 574 |
zab = zia - zib;
|
| 575 |
rab = sqrt(xab*xab + yab*yab + zab*zab);
|
| 576 |
xcb = xic - xib;
|
| 577 |
ycb = yic - yib;
|
| 578 |
zcb = zic - zib;
|
| 579 |
rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
|
| 580 |
xp = ycb*zab - zcb*yab;
|
| 581 |
yp = zcb*xab - xcb*zab;
|
| 582 |
zp = xcb*yab - ycb*xab;
|
| 583 |
rp = sqrt(xp*xp + yp*yp + zp*zp);
|
| 584 |
if (rp != 0.0)
|
| 585 |
{
|
| 586 |
dot = xab*xcb + yab*ycb + zab*zcb;
|
| 587 |
cosine = dot / (rab*rcb);
|
| 588 |
if (cosine < -1.0)
|
| 589 |
cosine = -1.0;
|
| 590 |
if (cosine > 1.0)
|
| 591 |
cosine = 1.0;
|
| 592 |
angle = radian * acos(cosine);
|
| 593 |
|
| 594 |
if (angles.angtype[i] == HARMONIC)
|
| 595 |
{
|
| 596 |
dt = angle - angles.anat[i];
|
| 597 |
dt2 = dt * dt;
|
| 598 |
dt3 = dt2 * dt;
|
| 599 |
dt4 = dt3 * dt;
|
| 600 |
deddt = units.angunit * angles.acon[i] * dt
|
| 601 |
* (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
|
| 602 |
+ 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
|
| 603 |
d2eddt2 = units.angunit * angles.acon[i]
|
| 604 |
* (2.0 + 6.0*units.cang*dt + 12.0*units.qang*dt2
|
| 605 |
+ 20.0*units.pang*dt3 + 30.0*units.sang*dt4);
|
| 606 |
terma = -radian/ (rab*rab*rp);
|
| 607 |
termc = radian / (rcb*rcb*rp);
|
| 608 |
} else
|
| 609 |
{
|
| 610 |
factor = angle/radian;
|
| 611 |
sine = sin(factor);
|
| 612 |
deddt = -units.angunit*angles.fcon[i]*
|
| 613 |
(angles.c1[i]*sine + 2*angles.c2[i]*sin(2*factor)
|
| 614 |
+ 3*angles.c3[i]*sin(3*factor) + 4*angles.c4[i]*sin(4*factor)
|
| 615 |
+ 5*angles.c5[i]*sin(5*factor) + 6*angles.c6[i]*sin(6*factor));
|
| 616 |
d2eddt2 = -units.angunit*angles.fcon[i]*
|
| 617 |
(angles.c1[i]*cosine + 4*angles.c2[i]*cos(2*factor)
|
| 618 |
+ 9*angles.c3[i]*cos(3*factor) + 16*angles.c4[i]*cos(4*factor)
|
| 619 |
+ 25*angles.c5[i]*cos(5*factor) + 36*angles.c6[i]*cos(6*factor));
|
| 620 |
terma = -1.0 / (rab*rab*rp);
|
| 621 |
termc = 1.0 / (rcb*rcb*rp);
|
| 622 |
}
|
| 623 |
ddtdxia = terma * (yab*zp-zab*yp);
|
| 624 |
ddtdyia = terma * (zab*xp-xab*zp);
|
| 625 |
ddtdzia = terma * (xab*yp-yab*xp);
|
| 626 |
ddtdxic = termc * (ycb*zp-zcb*yp);
|
| 627 |
ddtdyic = termc * (zcb*xp-xcb*zp);
|
| 628 |
ddtdzic = termc * (xcb*yp-ycb*xp);
|
| 629 |
ddtdxib = -ddtdxia - ddtdxic;
|
| 630 |
ddtdyib = -ddtdyia - ddtdyic;
|
| 631 |
ddtdzib = -ddtdzia - ddtdzic;
|
| 632 |
rab2 = 2.0 / (rab*rab);
|
| 633 |
xrab = xab * rab2;
|
| 634 |
yrab = yab * rab2;
|
| 635 |
zrab = zab * rab2;
|
| 636 |
rcb2 = 2.0 / (rcb*rcb);
|
| 637 |
xrcb = xcb * rcb2;
|
| 638 |
yrcb = ycb * rcb2;
|
| 639 |
zrcb = zcb * rcb2;
|
| 640 |
rp2 = 1.0 / (rp*rp);
|
| 641 |
xabp = (yab*zp-zab*yp) * rp2;
|
| 642 |
yabp = (zab*xp-xab*zp) * rp2;
|
| 643 |
zabp = (xab*yp-yab*xp) * rp2;
|
| 644 |
xcbp = (ycb*zp-zcb*yp) * rp2;
|
| 645 |
ycbp = (zcb*xp-xcb*zp) * rp2;
|
| 646 |
zcbp = (xcb*yp-ycb*xp) * rp2;
|
| 647 |
|
| 648 |
dxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab);
|
| 649 |
dxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab);
|
| 650 |
dxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab);
|
| 651 |
dyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab);
|
| 652 |
dyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab);
|
| 653 |
dziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab);
|
| 654 |
dxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb);
|
| 655 |
dxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb);
|
| 656 |
dxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb);
|
| 657 |
dyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb);
|
| 658 |
dyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb);
|
| 659 |
dziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb);
|
| 660 |
dxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp;
|
| 661 |
dxiayic = -terma*xab*yab - ddtdxia*yabp;
|
| 662 |
dxiazic = -terma*xab*zab - ddtdxia*zabp;
|
| 663 |
dyiaxic = -terma*xab*yab - ddtdyia*xabp;
|
| 664 |
dyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp;
|
| 665 |
dyiazic = -terma*yab*zab - ddtdyia*zabp;
|
| 666 |
dziaxic = -terma*xab*zab - ddtdzia*xabp;
|
| 667 |
dziayic = -terma*yab*zab - ddtdzia*yabp;
|
| 668 |
dziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp;
|
| 669 |
|
| 670 |
dxibxia = -dxiaxia - dxiaxic;
|
| 671 |
dxibyia = -dxiayia - dyiaxic;
|
| 672 |
dxibzia = -dxiazia - dziaxic;
|
| 673 |
dyibxia = -dxiayia - dxiayic;
|
| 674 |
dyibyia = -dyiayia - dyiayic;
|
| 675 |
dyibzia = -dyiazia - dziayic;
|
| 676 |
dzibxia = -dxiazia - dxiazic;
|
| 677 |
dzibyia = -dyiazia - dyiazic;
|
| 678 |
dzibzia = -dziazia - dziazic;
|
| 679 |
dxibxic = -dxicxic - dxiaxic;
|
| 680 |
dxibyic = -dxicyic - dxiayic;
|
| 681 |
dxibzic = -dxiczic - dxiazic;
|
| 682 |
dyibxic = -dxicyic - dyiaxic;
|
| 683 |
dyibyic = -dyicyic - dyiayic;
|
| 684 |
dyibzic = -dyiczic - dyiazic;
|
| 685 |
dzibxic = -dxiczic - dziaxic;
|
| 686 |
dzibyic = -dyiczic - dziayic;
|
| 687 |
dzibzic = -dziczic - dziazic;
|
| 688 |
dxibxib = -dxibxia - dxibxic;
|
| 689 |
dxibyib = -dxibyia - dxibyic;
|
| 690 |
dxibzib = -dxibzia - dxibzic;
|
| 691 |
dyibyib = -dyibyia - dyibyic;
|
| 692 |
dyibzib = -dyibzia - dyibzic;
|
| 693 |
dzibzib = -dzibzia - dzibzic;
|
| 694 |
|
| 695 |
if (ia == iatom)
|
| 696 |
{
|
| 697 |
hess.hessx[ia][0] += deddt*dxiaxia + d2eddt2*ddtdxia*ddtdxia;
|
| 698 |
hess.hessx[ia][1] += deddt*dxiayia + d2eddt2*ddtdxia*ddtdyia;
|
| 699 |
hess.hessx[ia][2] += deddt*dxiazia + d2eddt2*ddtdxia*ddtdzia;
|
| 700 |
hess.hessy[ia][0] += deddt*dxiayia + d2eddt2*ddtdyia*ddtdxia;
|
| 701 |
hess.hessy[ia][1] += deddt*dyiayia + d2eddt2*ddtdyia*ddtdyia;
|
| 702 |
hess.hessy[ia][2] += deddt*dyiazia + d2eddt2*ddtdyia*ddtdzia;
|
| 703 |
hess.hessz[ia][0] += deddt*dxiazia + d2eddt2*ddtdzia*ddtdxia;
|
| 704 |
hess.hessz[ia][1] += deddt*dyiazia + d2eddt2*ddtdzia*ddtdyia;
|
| 705 |
hess.hessz[ia][2] += deddt*dziazia + d2eddt2*ddtdzia*ddtdzia;
|
| 706 |
hess.hessx[ib][0] += deddt*dxibxia + d2eddt2*ddtdxia*ddtdxib;
|
| 707 |
hess.hessx[ib][1] += deddt*dyibxia + d2eddt2*ddtdxia*ddtdyib;
|
| 708 |
hess.hessx[ib][2] += deddt*dzibxia + d2eddt2*ddtdxia*ddtdzib;
|
| 709 |
hess.hessy[ib][0] += deddt*dxibyia + d2eddt2*ddtdyia*ddtdxib;
|
| 710 |
hess.hessy[ib][1] += deddt*dyibyia + d2eddt2*ddtdyia*ddtdyib;
|
| 711 |
hess.hessy[ib][2] += deddt*dzibyia + d2eddt2*ddtdyia*ddtdzib;
|
| 712 |
hess.hessz[ib][0] += deddt*dxibzia + d2eddt2*ddtdzia*ddtdxib;
|
| 713 |
hess.hessz[ib][1] += deddt*dyibzia + d2eddt2*ddtdzia*ddtdyib;
|
| 714 |
hess.hessz[ib][2] += deddt*dzibzia + d2eddt2*ddtdzia*ddtdzib;
|
| 715 |
hess.hessx[ic][0] += deddt*dxiaxic + d2eddt2*ddtdxia*ddtdxic;
|
| 716 |
hess.hessx[ic][1] += deddt*dxiayic + d2eddt2*ddtdxia*ddtdyic;
|
| 717 |
hess.hessx[ic][2] += deddt*dxiazic + d2eddt2*ddtdxia*ddtdzic;
|
| 718 |
hess.hessy[ic][0] += deddt*dyiaxic + d2eddt2*ddtdyia*ddtdxic;
|
| 719 |
hess.hessy[ic][1] += deddt*dyiayic + d2eddt2*ddtdyia*ddtdyic;
|
| 720 |
hess.hessy[ic][2] += deddt*dyiazic + d2eddt2*ddtdyia*ddtdzic;
|
| 721 |
hess.hessz[ic][0] += deddt*dziaxic + d2eddt2*ddtdzia*ddtdxic;
|
| 722 |
hess.hessz[ic][1] += deddt*dziayic + d2eddt2*ddtdzia*ddtdyic;
|
| 723 |
hess.hessz[ic][2] += deddt*dziazic + d2eddt2*ddtdzia*ddtdzic;
|
| 724 |
} else if (ib == iatom)
|
| 725 |
{
|
| 726 |
hess.hessx[ib][0] += deddt*dxibxib + d2eddt2*ddtdxib*ddtdxib;
|
| 727 |
hess.hessx[ib][1] += deddt*dxibyib + d2eddt2*ddtdxib*ddtdyib;
|
| 728 |
hess.hessx[ib][2] += deddt*dxibzib + d2eddt2*ddtdxib*ddtdzib;
|
| 729 |
hess.hessy[ib][0] += deddt*dxibyib + d2eddt2*ddtdyib*ddtdxib;
|
| 730 |
hess.hessy[ib][1] += deddt*dyibyib + d2eddt2*ddtdyib*ddtdyib;
|
| 731 |
hess.hessy[ib][2] += deddt*dyibzib + d2eddt2*ddtdyib*ddtdzib;
|
| 732 |
hess.hessz[ib][0] += deddt*dxibzib + d2eddt2*ddtdzib*ddtdxib;
|
| 733 |
hess.hessz[ib][1] += deddt*dyibzib + d2eddt2*ddtdzib*ddtdyib;
|
| 734 |
hess.hessz[ib][2] += deddt*dzibzib + d2eddt2*ddtdzib*ddtdzib;
|
| 735 |
hess.hessx[ia][0] += deddt*dxibxia + d2eddt2*ddtdxib*ddtdxia;
|
| 736 |
hess.hessx[ia][1] += deddt*dxibyia + d2eddt2*ddtdxib*ddtdyia;
|
| 737 |
hess.hessx[ia][2] += deddt*dxibzia + d2eddt2*ddtdxib*ddtdzia;
|
| 738 |
hess.hessy[ia][0] += deddt*dyibxia + d2eddt2*ddtdyib*ddtdxia;
|
| 739 |
hess.hessy[ia][1] += deddt*dyibyia + d2eddt2*ddtdyib*ddtdyia;
|
| 740 |
hess.hessy[ia][2] += deddt*dyibzia + d2eddt2*ddtdyib*ddtdzia;
|
| 741 |
hess.hessz[ia][0] += deddt*dzibxia + d2eddt2*ddtdzib*ddtdxia;
|
| 742 |
hess.hessz[ia][1] += deddt*dzibyia + d2eddt2*ddtdzib*ddtdyia;
|
| 743 |
hess.hessz[ia][2] += deddt*dzibzia + d2eddt2*ddtdzib*ddtdzia;
|
| 744 |
hess.hessx[ic][0] += deddt*dxibxic + d2eddt2*ddtdxib*ddtdxic;
|
| 745 |
hess.hessx[ic][1] += deddt*dxibyic + d2eddt2*ddtdxib*ddtdyic;
|
| 746 |
hess.hessx[ic][2] += deddt*dxibzic + d2eddt2*ddtdxib*ddtdzic;
|
| 747 |
hess.hessy[ic][0] += deddt*dyibxic + d2eddt2*ddtdyib*ddtdxic;
|
| 748 |
hess.hessy[ic][1] += deddt*dyibyic + d2eddt2*ddtdyib*ddtdyic;
|
| 749 |
hess.hessy[ic][2] += deddt*dyibzic + d2eddt2*ddtdyib*ddtdzic;
|
| 750 |
hess.hessz[ic][0] += deddt*dzibxic + d2eddt2*ddtdzib*ddtdxic;
|
| 751 |
hess.hessz[ic][1] += deddt*dzibyic + d2eddt2*ddtdzib*ddtdyic;
|
| 752 |
hess.hessz[ic][2] += deddt*dzibzic + d2eddt2*ddtdzib*ddtdzic;
|
| 753 |
}else if (ic == iatom)
|
| 754 |
{
|
| 755 |
hess.hessx[ic][0] += deddt*dxicxic + d2eddt2*ddtdxic*ddtdxic;
|
| 756 |
hess.hessx[ic][1] += deddt*dxicyic + d2eddt2*ddtdxic*ddtdyic;
|
| 757 |
hess.hessx[ic][2] += deddt*dxiczic + d2eddt2*ddtdxic*ddtdzic;
|
| 758 |
hess.hessy[ic][0] += deddt*dxicyic + d2eddt2*ddtdyic*ddtdxic;
|
| 759 |
hess.hessy[ic][1] += deddt*dyicyic + d2eddt2*ddtdyic*ddtdyic;
|
| 760 |
hess.hessy[ic][2] += deddt*dyiczic + d2eddt2*ddtdyic*ddtdzic;
|
| 761 |
hess.hessz[ic][0] += deddt*dxiczic + d2eddt2*ddtdzic*ddtdxic;
|
| 762 |
hess.hessz[ic][1] += deddt*dyiczic + d2eddt2*ddtdzic*ddtdyic;
|
| 763 |
hess.hessz[ic][2] += deddt*dziczic + d2eddt2*ddtdzic*ddtdzic;
|
| 764 |
hess.hessx[ib][0] += deddt*dxibxic + d2eddt2*ddtdxic*ddtdxib;
|
| 765 |
hess.hessx[ib][1] += deddt*dyibxic + d2eddt2*ddtdxic*ddtdyib;
|
| 766 |
hess.hessx[ib][2] += deddt*dzibxic + d2eddt2*ddtdxic*ddtdzib;
|
| 767 |
hess.hessy[ib][0] += deddt*dxibyic + d2eddt2*ddtdyic*ddtdxib;
|
| 768 |
hess.hessy[ib][1] += deddt*dyibyic + d2eddt2*ddtdyic*ddtdyib;
|
| 769 |
hess.hessy[ib][2] += deddt*dzibyic + d2eddt2*ddtdyic*ddtdzib;
|
| 770 |
hess.hessz[ib][0] += deddt*dxibzic + d2eddt2*ddtdzic*ddtdxib;
|
| 771 |
hess.hessz[ib][1] += deddt*dyibzic + d2eddt2*ddtdzic*ddtdyib;
|
| 772 |
hess.hessz[ib][2] += deddt*dzibzic + d2eddt2*ddtdzic*ddtdzib;
|
| 773 |
hess.hessx[ia][0] += deddt*dxiaxic + d2eddt2*ddtdxic*ddtdxia;
|
| 774 |
hess.hessx[ia][1] += deddt*dyiaxic + d2eddt2*ddtdxic*ddtdyia;
|
| 775 |
hess.hessx[ia][2] += deddt*dziaxic + d2eddt2*ddtdxic*ddtdzia;
|
| 776 |
hess.hessy[ia][0] += deddt*dxiayic + d2eddt2*ddtdyic*ddtdxia;
|
| 777 |
hess.hessy[ia][1] += deddt*dyiayic + d2eddt2*ddtdyic*ddtdyia;
|
| 778 |
hess.hessy[ia][2] += deddt*dziayic + d2eddt2*ddtdyic*ddtdzia;
|
| 779 |
hess.hessz[ia][0] += deddt*dxiazic + d2eddt2*ddtdzic*ddtdxia;
|
| 780 |
hess.hessz[ia][1] += deddt*dyiazic + d2eddt2*ddtdzic*ddtdyia;
|
| 781 |
hess.hessz[ia][2] += deddt*dziazic + d2eddt2*ddtdzic*ddtdzia;
|
| 782 |
}
|
| 783 |
}
|
| 784 |
}
|
| 785 |
}
|
| 786 |
}
|
| 787 |
}
|
| 788 |
// ================================================================
|
| 789 |
void eangle2b(int i)
|
| 790 |
{
|
| 791 |
int ia,ib,ic,id;
|
| 792 |
double angle,dot,cosine;
|
| 793 |
double dt,dt2,dt3,dt4;
|
| 794 |
double deddt,terma,termc,term;
|
| 795 |
double xia,yia,zia;
|
| 796 |
double xib,yib,zib;
|
| 797 |
double xic,yic,zic;
|
| 798 |
double xid,yid,zid;
|
| 799 |
double xad,yad,zad;
|
| 800 |
double xbd,ybd,zbd;
|
| 801 |
double xcd,ycd,zcd;
|
| 802 |
double xip,yip,zip;
|
| 803 |
double xap,yap,zap,rap2;
|
| 804 |
double xcp,ycp,zcp,rcp2;
|
| 805 |
double xt,yt,zt,rt2,ptrt2;
|
| 806 |
double xm,ym,zm,rm,delta,delta2;
|
| 807 |
double dedxia,dedyia,dedzia;
|
| 808 |
double dedxib,dedyib,dedzib;
|
| 809 |
double dedxic,dedyic,dedzic;
|
| 810 |
double dedxid,dedyid,dedzid;
|
| 811 |
double dedxip,dedyip,dedzip;
|
| 812 |
double dpdxia,dpdyia,dpdzia;
|
| 813 |
double dpdxic,dpdyic,dpdzic;
|
| 814 |
|
| 815 |
ia = angles.i13[i][0];
|
| 816 |
ib = angles.i13[i][1];
|
| 817 |
ic = angles.i13[i][2];
|
| 818 |
id = angles.i13[i][3];
|
| 819 |
xia = atom[ia].x;
|
| 820 |
yia = atom[ia].y;
|
| 821 |
zia = atom[ia].z;
|
| 822 |
xib = atom[ib].x;
|
| 823 |
yib = atom[ib].y;
|
| 824 |
zib = atom[ib].z;
|
| 825 |
xic = atom[ic].x;
|
| 826 |
yic = atom[ic].y;
|
| 827 |
zic = atom[ic].z;
|
| 828 |
xid = atom[id].x;
|
| 829 |
yid = atom[id].y;
|
| 830 |
zid = atom[id].z;
|
| 831 |
|
| 832 |
deriv.dea[ia][0] = 0.0;
|
| 833 |
deriv.dea[ia][1] = 0.0;
|
| 834 |
deriv.dea[ia][2] = 0.0;
|
| 835 |
deriv.dea[ib][0] = 0.0;
|
| 836 |
deriv.dea[ib][1] = 0.0;
|
| 837 |
deriv.dea[ib][2] = 0.0;
|
| 838 |
deriv.dea[ic][0] = 0.0;
|
| 839 |
deriv.dea[ic][1] = 0.0;
|
| 840 |
deriv.dea[ic][2] = 0.0;
|
| 841 |
deriv.dea[id][0] = 0.0;
|
| 842 |
deriv.dea[id][1] = 0.0;
|
| 843 |
deriv.dea[id][2] = 0.0;
|
| 844 |
|
| 845 |
xad = xia - xid;
|
| 846 |
yad = yia - yid;
|
| 847 |
zad = zia - zid;
|
| 848 |
xbd = xib - xid;
|
| 849 |
ybd = yib - yid;
|
| 850 |
zbd = zib - zid;
|
| 851 |
xcd = xic - xid;
|
| 852 |
ycd = yic - yid;
|
| 853 |
zcd = zic - zid;
|
| 854 |
xt = yad*zcd - zad*ycd;
|
| 855 |
yt = zad*xcd - xad*zcd;
|
| 856 |
zt = xad*ycd - yad*xcd;
|
| 857 |
rt2 = xt*xt + yt*yt + zt*zt;
|
| 858 |
delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
|
| 859 |
xip = xib + xt*delta;
|
| 860 |
yip = yib + yt*delta;
|
| 861 |
zip = zib + zt*delta;
|
| 862 |
xap = xia - xip;
|
| 863 |
yap = yia - yip;
|
| 864 |
zap = zia - zip;
|
| 865 |
rap2 = xap*xap + yap*yap + zap*zap;
|
| 866 |
xcp = xic - xip;
|
| 867 |
ycp = yic - yip;
|
| 868 |
zcp = zic - zip;
|
| 869 |
rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
|
| 870 |
xm = ycp*zap - zcp*yap;
|
| 871 |
ym = zcp*xap - xcp*zap;
|
| 872 |
zm = xcp*yap - ycp*xap;
|
| 873 |
rm = sqrt(xm*xm + ym*ym + zm*zm);
|
| 874 |
if (rm != 0.0)
|
| 875 |
{
|
| 876 |
dot = xap*xcp + yap*ycp + zap*zcp;
|
| 877 |
cosine = dot / sqrt(rap2*rcp2);
|
| 878 |
if (cosine < -1.0)
|
| 879 |
cosine = -1.0;
|
| 880 |
if (cosine > 1.0)
|
| 881 |
cosine = 1.0;
|
| 882 |
angle = radian * acos(cosine);
|
| 883 |
|
| 884 |
dt = angle - angles.anat[i];
|
| 885 |
dt2 = dt * dt;
|
| 886 |
dt3 = dt2 * dt;
|
| 887 |
dt4 = dt2 * dt2;
|
| 888 |
deddt = units.angunit * angles.acon[i] * dt
|
| 889 |
* (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
|
| 890 |
+ 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
|
| 891 |
|
| 892 |
deddt = deddt * radian;
|
| 893 |
terma = -deddt / (rap2*rm);
|
| 894 |
termc = deddt / (rcp2*rm);
|
| 895 |
dedxia = terma * (yap*zm-zap*ym);
|
| 896 |
dedyia = terma * (zap*xm-xap*zm);
|
| 897 |
dedzia = terma * (xap*ym-yap*xm);
|
| 898 |
dedxic = termc * (ycp*zm-zcp*ym);
|
| 899 |
dedyic = termc * (zcp*xm-xcp*zm);
|
| 900 |
dedzic = termc * (xcp*ym-ycp*xm);
|
| 901 |
dedxip = -dedxia - dedxic;
|
| 902 |
dedyip = -dedyia - dedyic;
|
| 903 |
dedzip = -dedzia - dedzic;
|
| 904 |
|
| 905 |
delta2 = 2.0 * delta;
|
| 906 |
ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
|
| 907 |
term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
|
| 908 |
dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
|
| 909 |
term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
|
| 910 |
dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
|
| 911 |
term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
|
| 912 |
dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
|
| 913 |
term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
|
| 914 |
dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
|
| 915 |
term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
|
| 916 |
dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
|
| 917 |
term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
|
| 918 |
dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
|
| 919 |
dedxia = dedxia + dpdxia;
|
| 920 |
dedyia = dedyia + dpdyia;
|
| 921 |
dedzia = dedzia + dpdzia;
|
| 922 |
dedxib = dedxip;
|
| 923 |
dedyib = dedyip;
|
| 924 |
dedzib = dedzip;
|
| 925 |
dedxic = dedxic + dpdxic;
|
| 926 |
dedyic = dedyic + dpdyic;
|
| 927 |
dedzic = dedzic + dpdzic;
|
| 928 |
dedxid = -dedxia - dedxib - dedxic;
|
| 929 |
dedyid = -dedyia - dedyib - dedyic;
|
| 930 |
dedzid = -dedzia - dedzib - dedzic;
|
| 931 |
|
| 932 |
deriv.dea[ia][0] = dedxia;
|
| 933 |
deriv.dea[ia][1] = dedyia;
|
| 934 |
deriv.dea[ia][2] = dedzia;
|
| 935 |
deriv.dea[ib][0] = dedxib;
|
| 936 |
deriv.dea[ib][1] = dedyib;
|
| 937 |
deriv.dea[ib][2] = dedzib;
|
| 938 |
deriv.dea[ic][0] = dedxic;
|
| 939 |
deriv.dea[ic][1] = dedyic;
|
| 940 |
deriv.dea[ic][2] = dedzic;
|
| 941 |
deriv.dea[id][0] = dedxid;
|
| 942 |
deriv.dea[id][1] = dedyid;
|
| 943 |
deriv.dea[id][2] = dedzid;
|
| 944 |
}
|
| 945 |
}
|
| 946 |
|
