| 1 |
#define EXTERN extern |
| 2 |
|
| 3 |
#include "pcwin.h" |
| 4 |
#include "pcmod.h" |
| 5 |
|
| 6 |
#include "energies.h" |
| 7 |
#include "angles.h" |
| 8 |
#include "bonds_ff.h" |
| 9 |
#include "hess.h" |
| 10 |
#include "derivs.h" |
| 11 |
|
| 12 |
EXTERN struct t_strbnd { |
| 13 |
int nstrbnd, isb[MAXANG][3]; |
| 14 |
float ksb1[MAXANG],ksb2[MAXANG]; |
| 15 |
} strbnd; |
| 16 |
|
| 17 |
EXTERN struct t_minim_values { |
| 18 |
int iprint, ndc, nconst; |
| 19 |
float dielc; |
| 20 |
} minim_values; |
| 21 |
|
| 22 |
void estrbnd() |
| 23 |
{ |
| 24 |
int i,ij,j,k,ia,ib,ic; |
| 25 |
double e,dr,dt,angle,dot,cosine,dr1; |
| 26 |
double xia,yia,zia,xib,yib,zib; |
| 27 |
double xic,yic,zic; |
| 28 |
double rab,xab,yab,zab; |
| 29 |
double rcb,xcb,ycb,zcb; |
| 30 |
|
| 31 |
energies.estrbnd = 0; |
| 32 |
if (minim_values.iprint && strbnd.nstrbnd > 0) |
| 33 |
{ |
| 34 |
fprintf(pcmlogfile,"\nStretch-Bend Terms\n"); |
| 35 |
fprintf(pcmlogfile," At1 At2 At3 Angle Anat StbnCst1 StbnCst2 Estb\n"); |
| 36 |
} |
| 37 |
|
| 38 |
for (i=0; i < strbnd.nstrbnd; i++) |
| 39 |
{ |
| 40 |
ij = strbnd.isb[i][0]; |
| 41 |
ia = angles.i13[ij][0]; |
| 42 |
ib = angles.i13[ij][1]; |
| 43 |
ic = angles.i13[ij][2]; |
| 44 |
if (atom[ia].use || atom[ib].use || atom[ic].use) |
| 45 |
{ |
| 46 |
xia = atom[ia].x; |
| 47 |
yia = atom[ia].y; |
| 48 |
zia = atom[ia].z; |
| 49 |
xib = atom[ib].x; |
| 50 |
yib = atom[ib].y; |
| 51 |
zib = atom[ib].z; |
| 52 |
xic = atom[ic].x; |
| 53 |
yic = atom[ic].y; |
| 54 |
zic = atom[ic].z; |
| 55 |
|
| 56 |
xab = xia - xib; |
| 57 |
yab = yia - yib; |
| 58 |
zab = zia - zib; |
| 59 |
rab = sqrt(xab*xab + yab*yab + zab*zab); |
| 60 |
xcb = xic - xib; |
| 61 |
ycb = yic - yib; |
| 62 |
zcb = zic - zib; |
| 63 |
rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb); |
| 64 |
if (rab*rcb != 0.0) |
| 65 |
{ |
| 66 |
dot = xab*xcb + yab*ycb + zab*zcb; |
| 67 |
cosine = dot / (rab*rcb); |
| 68 |
if (cosine < -1.0) |
| 69 |
cosine = -1.0; |
| 70 |
if (cosine > 1.0) |
| 71 |
cosine = 1.0; |
| 72 |
angle = radian * acos(cosine); |
| 73 |
dt = angle - angles.anat[ij]; |
| 74 |
j = strbnd.isb[i][1]; |
| 75 |
k = strbnd.isb[i][2]; |
| 76 |
dr = 0.0; |
| 77 |
dr1 = 0.0; |
| 78 |
if (j > -1) dr = (rab - bonds_ff.bl[j])*strbnd.ksb1[i]; |
| 79 |
if (k > -1) dr1 = (rcb - bonds_ff.bl[k])*strbnd.ksb2[i]; |
| 80 |
e = units.stbnunit*dt*(dr+dr1); |
| 81 |
energies.estrbnd += e; |
| 82 |
atom[ia].energy += e; |
| 83 |
atom[ib].energy += e; |
| 84 |
atom[ic].energy += e; |
| 85 |
|
| 86 |
if (minim_values.iprint) |
| 87 |
fprintf(pcmlogfile,"StrBnd: %2s(%-3d)-%2s(%-3d)-%2s(%-3d) : %-8.3f %-8.3f %-8.3f %-8.3f = %-8.4f\n", |
| 88 |
atom[ia].name, ia,atom[ib].name, ib,atom[ic].name, ic,angle,angles.anat[ij], |
| 89 |
strbnd.ksb1[i],strbnd.ksb2[i], e); |
| 90 |
} |
| 91 |
} |
| 92 |
} |
| 93 |
} |
| 94 |
|
| 95 |
void estrbnd1() |
| 96 |
{ |
| 97 |
int i,ij,j,k,ia,ib,ic; |
| 98 |
double e,dr,dr1,dt,angle,dot,cosine; |
| 99 |
double xia,yia,zia,xib,yib,zib; |
| 100 |
double xic,yic,zic; |
| 101 |
double rab,xab,yab,zab; |
| 102 |
double rcb,xcb,ycb,zcb; |
| 103 |
double xp,yp,zp,rp; |
| 104 |
double dedxia,dedyia,dedzia; |
| 105 |
double dedxib,dedyib,dedzib; |
| 106 |
double dedxic,dedyic,dedzic; |
| 107 |
double term,terma,termc, rab2, rcb2; |
| 108 |
double ksb1, ksb2; |
| 109 |
double dtemp; |
| 110 |
|
| 111 |
energies.estrbnd = 0; |
| 112 |
for (i=0; i <= natom; i++) |
| 113 |
{ |
| 114 |
deriv.destb[i][0] = 0.0; |
| 115 |
deriv.destb[i][1] = 0.0; |
| 116 |
deriv.destb[i][2] = 0.0; |
| 117 |
} |
| 118 |
|
| 119 |
for (i=0; i < strbnd.nstrbnd; i++) |
| 120 |
{ |
| 121 |
ij = strbnd.isb[i][0]; |
| 122 |
ia = angles.i13[ij][0]; |
| 123 |
ib = angles.i13[ij][1]; |
| 124 |
ic = angles.i13[ij][2]; |
| 125 |
if (atom[ia].use || atom[ib].use || atom[ic].use) |
| 126 |
{ |
| 127 |
xia = atom[ia].x; |
| 128 |
yia = atom[ia].y; |
| 129 |
zia = atom[ia].z; |
| 130 |
xib = atom[ib].x; |
| 131 |
yib = atom[ib].y; |
| 132 |
zib = atom[ib].z; |
| 133 |
xic = atom[ic].x; |
| 134 |
yic = atom[ic].y; |
| 135 |
zic = atom[ic].z; |
| 136 |
|
| 137 |
xab = xia - xib; |
| 138 |
yab = yia - yib; |
| 139 |
zab = zia - zib; |
| 140 |
rab2 = xab*xab + yab*yab + zab*zab; |
| 141 |
rab = sqrt(xab*xab + yab*yab + zab*zab); |
| 142 |
xcb = xic - xib; |
| 143 |
ycb = yic - yib; |
| 144 |
zcb = zic - zib; |
| 145 |
rcb2 = xcb*xcb + ycb*ycb + zcb*zcb; |
| 146 |
rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb); |
| 147 |
xp = ycb*zab - zcb*yab; |
| 148 |
yp = zcb*xab - xcb*zab; |
| 149 |
zp = xcb*yab - ycb*xab; |
| 150 |
rp = sqrt(xp*xp + yp*yp + zp*zp); |
| 151 |
if (rp != 0.0) |
| 152 |
{ |
| 153 |
dot = xab*xcb + yab*ycb + zab*zcb; |
| 154 |
cosine = dot / (rab*rcb); |
| 155 |
if (cosine < -1.0) |
| 156 |
cosine = -1.0; |
| 157 |
if (cosine > 1.0) |
| 158 |
cosine = 1.0; |
| 159 |
angle = radian * acos(cosine); |
| 160 |
dt = angle - angles.anat[ij]; |
| 161 |
|
| 162 |
dtemp = dot*dot/(rab2*rcb2); |
| 163 |
if (fabs(dtemp) >= 1.0) |
| 164 |
dtemp = 0.9999*dtemp/fabs(dtemp); |
| 165 |
termc = sqrt(1.00 - dtemp); |
| 166 |
terma = rab*rcb; |
| 167 |
j = strbnd.isb[i][1]; |
| 168 |
k = strbnd.isb[i][2]; |
| 169 |
term = -units.stbnunit*radian; |
| 170 |
if (j < 0) |
| 171 |
{ |
| 172 |
dr = 0.0; |
| 173 |
ksb1 = 0.0; |
| 174 |
} else |
| 175 |
{ |
| 176 |
dr = strbnd.ksb1[i]*(rab - bonds_ff.bl[j]); |
| 177 |
ksb1 = strbnd.ksb1[i]; |
| 178 |
} |
| 179 |
if (k < 0) |
| 180 |
{ |
| 181 |
dr1 = 0.0; |
| 182 |
ksb2 = 0.0; |
| 183 |
} else |
| 184 |
{ |
| 185 |
dr1 = strbnd.ksb2[i]*(rcb - bonds_ff.bl[k]); |
| 186 |
ksb2 = strbnd.ksb2[i]; |
| 187 |
} |
| 188 |
e = units.stbnunit*dt*(dr+dr1); |
| 189 |
dedxia = (term *(xcb/terma - (xab*dot)/(rab2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb1*xab)/rab; |
| 190 |
dedyia = (term *(ycb/terma - (yab*dot)/(rab2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb1*yab)/rab; |
| 191 |
dedzia = (term *(zcb/terma - (zab*dot)/(rab2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb1*zab)/rab; |
| 192 |
|
| 193 |
dedxic = (term *(xab/terma - (xcb*dot)/(rcb2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb2*xcb)/rcb; |
| 194 |
dedyic = (term *(yab/terma - (ycb*dot)/(rcb2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb2*ycb)/rcb; |
| 195 |
dedzic = (term *(zab/terma - (zcb*dot)/(rcb2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb2*zcb)/rcb; |
| 196 |
|
| 197 |
dedxib = ( (term /(rab*rcb))*( (xcb*dot)/rcb2 + (xab*dot)/rab2 + (2.0*xib-xia-xic))*(dr+dr1))/termc + |
| 198 |
units.stbnunit*dt*( -(ksb1*xab)/rab - ksb2*xcb/rcb); |
| 199 |
dedyib = ( (term /(rab*rcb))*( (ycb*dot)/rcb2 + (yab*dot)/rab2 + (2.0*yib-yia-yic))*(dr+dr1))/termc + |
| 200 |
units.stbnunit*dt*( -(ksb1*yab)/rab - ksb2*ycb/rcb); |
| 201 |
dedzib = ( (term /(rab*rcb))*( (zcb*dot)/rcb2 + (zab*dot)/rab2 + (2.0*zib-zia-zic))*(dr+dr1))/termc + |
| 202 |
units.stbnunit*dt*( -(ksb1*zab)/rab - ksb2*zcb/rcb); |
| 203 |
|
| 204 |
energies.estrbnd += e; |
| 205 |
deriv.destb[ia][0] += dedxia; |
| 206 |
deriv.destb[ia][1] += dedyia; |
| 207 |
deriv.destb[ia][2] += dedzia; |
| 208 |
|
| 209 |
deriv.destb[ib][0] += dedxib; |
| 210 |
deriv.destb[ib][1] += dedyib; |
| 211 |
deriv.destb[ib][2] += dedzib; |
| 212 |
|
| 213 |
deriv.destb[ic][0] += dedxic; |
| 214 |
deriv.destb[ic][1] += dedyic; |
| 215 |
deriv.destb[ic][2] += dedzic; |
| 216 |
|
| 217 |
} |
| 218 |
} |
| 219 |
} |
| 220 |
} |
| 221 |
|
| 222 |
void estrbnd2(int iatom) |
| 223 |
{ |
| 224 |
int i,j,k,ij; |
| 225 |
int ia,ib,ic; |
| 226 |
double dt,dr,angle,dot,cosine; |
| 227 |
double xia,yia,zia; |
| 228 |
double xib,yib,zib; |
| 229 |
double xic,yic,zic; |
| 230 |
double xab,yab,zab,rab; |
| 231 |
double xcb,ycb,zcb,rcb; |
| 232 |
double xp,yp,zp,rp,rp2; |
| 233 |
double terma,termc; |
| 234 |
double xrab,yrab,zrab,rab2; |
| 235 |
double xrcb,yrcb,zrcb,rcb2; |
| 236 |
double xabp,yabp,zabp; |
| 237 |
double xcbp,ycbp,zcbp; |
| 238 |
double ddtdxia,ddtdyia,ddtdzia; |
| 239 |
double ddtdxib,ddtdyib,ddtdzib; |
| 240 |
double ddtdxic,ddtdyic,ddtdzic; |
| 241 |
double ddrdxia,ddrdyia,ddrdzia; |
| 242 |
double ddrdxib,ddrdyib,ddrdzib; |
| 243 |
double ddrdxic,ddrdyic,ddrdzic; |
| 244 |
double dtxiaxia,dtxiayia,dtxiazia; |
| 245 |
double dtxibxib,dtxibyib,dtxibzib; |
| 246 |
double dtxicxic,dtxicyic,dtxiczic; |
| 247 |
double dtyiayia,dtyiazia,dtziazia; |
| 248 |
double dtyibyib,dtyibzib,dtzibzib; |
| 249 |
double dtyicyic,dtyiczic,dtziczic; |
| 250 |
double dtxibxia,dtxibyia,dtxibzia; |
| 251 |
double dtyibxia,dtyibyia,dtyibzia; |
| 252 |
double dtzibxia,dtzibyia,dtzibzia; |
| 253 |
double dtxibxic,dtxibyic,dtxibzic; |
| 254 |
double dtyibxic,dtyibyic,dtyibzic; |
| 255 |
double dtzibxic,dtzibyic,dtzibzic; |
| 256 |
double dtxiaxic,dtxiayic,dtxiazic; |
| 257 |
double dtyiaxic,dtyiayic,dtyiazic; |
| 258 |
double dtziaxic,dtziayic,dtziazic; |
| 259 |
double drxiaxia,drxiayia,drxiazia; |
| 260 |
double drxibxib,drxibyib,drxibzib; |
| 261 |
double drxicxic,drxicyic,drxiczic; |
| 262 |
double dryiayia,dryiazia,drziazia; |
| 263 |
double dryibyib,dryibzib,drzibzib; |
| 264 |
double dryicyic,dryiczic,drziczic; |
| 265 |
|
| 266 |
for (i=0; i < strbnd.nstrbnd; i++) |
| 267 |
{ |
| 268 |
ij = strbnd.isb[i][0]; |
| 269 |
ia = angles.i13[ij][0]; |
| 270 |
ib = angles.i13[ij][1]; |
| 271 |
ic = angles.i13[ij][2]; |
| 272 |
if (iatom == ia || iatom == ib || iatom == ic) |
| 273 |
{ |
| 274 |
xia = atom[ia].x; |
| 275 |
yia = atom[ia].y; |
| 276 |
zia = atom[ia].z; |
| 277 |
xib = atom[ib].x; |
| 278 |
yib = atom[ib].y; |
| 279 |
zib = atom[ib].z; |
| 280 |
xic = atom[ic].x; |
| 281 |
yic = atom[ic].y; |
| 282 |
zic = atom[ic].z; |
| 283 |
xab = xia - xib; |
| 284 |
yab = yia - yib; |
| 285 |
zab = zia - zib; |
| 286 |
rab = sqrt(xab*xab + yab*yab + zab*zab); |
| 287 |
xcb = xic - xib; |
| 288 |
ycb = yic - yib; |
| 289 |
zcb = zic - zib; |
| 290 |
rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb); |
| 291 |
xp = ycb*zab - zcb*yab; |
| 292 |
yp = zcb*xab - xcb*zab; |
| 293 |
zp = xcb*yab - ycb*xab; |
| 294 |
rp = sqrt(xp*xp + yp*yp + zp*zp); |
| 295 |
if (rp != 0.0) |
| 296 |
{ |
| 297 |
dot = xab*xcb + yab*ycb + zab*zcb; |
| 298 |
cosine = dot / (rab*rcb); |
| 299 |
if (cosine < -1.0) |
| 300 |
cosine = -1.0; |
| 301 |
if (cosine > 1.0) |
| 302 |
cosine = 1.0; |
| 303 |
angle = radian * acos(cosine); |
| 304 |
|
| 305 |
dt = angle - angles.anat[ij]; |
| 306 |
terma = -radian / (rab*rab*rp); |
| 307 |
termc = radian / (rcb*rcb*rp); |
| 308 |
ddtdxia = terma * (yab*zp-zab*yp); |
| 309 |
ddtdyia = terma * (zab*xp-xab*zp); |
| 310 |
ddtdzia = terma * (xab*yp-yab*xp); |
| 311 |
ddtdxic = termc * (ycb*zp-zcb*yp); |
| 312 |
ddtdyic = termc * (zcb*xp-xcb*zp); |
| 313 |
ddtdzic = termc * (xcb*yp-ycb*xp); |
| 314 |
ddtdxib = -ddtdxia - ddtdxic; |
| 315 |
ddtdyib = -ddtdyia - ddtdyic; |
| 316 |
ddtdzib = -ddtdzia - ddtdzic; |
| 317 |
|
| 318 |
rab2 = 2.0 / (rab*rab); |
| 319 |
xrab = xab * rab2; |
| 320 |
yrab = yab * rab2; |
| 321 |
zrab = zab * rab2; |
| 322 |
rcb2 = 2.0 / (rcb*rcb); |
| 323 |
xrcb = xcb * rcb2; |
| 324 |
yrcb = ycb * rcb2; |
| 325 |
zrcb = zcb * rcb2; |
| 326 |
rp2 = 1.0 / (rp*rp); |
| 327 |
xabp = (yab*zp-zab*yp) * rp2; |
| 328 |
yabp = (zab*xp-xab*zp) * rp2; |
| 329 |
zabp = (xab*yp-yab*xp) * rp2; |
| 330 |
xcbp = (ycb*zp-zcb*yp) * rp2; |
| 331 |
ycbp = (zcb*xp-xcb*zp) * rp2; |
| 332 |
zcbp = (xcb*yp-ycb*xp) * rp2; |
| 333 |
|
| 334 |
dtxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab); |
| 335 |
dtxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab); |
| 336 |
dtxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab); |
| 337 |
dtyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab); |
| 338 |
dtyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab); |
| 339 |
dtziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab); |
| 340 |
dtxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb); |
| 341 |
dtxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb); |
| 342 |
dtxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb); |
| 343 |
dtyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb); |
| 344 |
dtyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb); |
| 345 |
dtziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb); |
| 346 |
dtxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp; |
| 347 |
dtxiayic = -terma*xab*yab - ddtdxia*yabp; |
| 348 |
dtxiazic = -terma*xab*zab - ddtdxia*zabp; |
| 349 |
dtyiaxic = -terma*xab*yab - ddtdyia*xabp; |
| 350 |
dtyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp; |
| 351 |
dtyiazic = -terma*yab*zab - ddtdyia*zabp; |
| 352 |
dtziaxic = -terma*xab*zab - ddtdzia*xabp; |
| 353 |
dtziayic = -terma*yab*zab - ddtdzia*yabp; |
| 354 |
dtziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp; |
| 355 |
|
| 356 |
dtxibxia = -dtxiaxia - dtxiaxic; |
| 357 |
dtxibyia = -dtxiayia - dtyiaxic; |
| 358 |
dtxibzia = -dtxiazia - dtziaxic; |
| 359 |
dtyibxia = -dtxiayia - dtxiayic; |
| 360 |
dtyibyia = -dtyiayia - dtyiayic; |
| 361 |
dtyibzia = -dtyiazia - dtziayic; |
| 362 |
dtzibxia = -dtxiazia - dtxiazic; |
| 363 |
dtzibyia = -dtyiazia - dtyiazic; |
| 364 |
dtzibzia = -dtziazia - dtziazic; |
| 365 |
dtxibxic = -dtxicxic - dtxiaxic; |
| 366 |
dtxibyic = -dtxicyic - dtxiayic; |
| 367 |
dtxibzic = -dtxiczic - dtxiazic; |
| 368 |
dtyibxic = -dtxicyic - dtyiaxic; |
| 369 |
dtyibyic = -dtyicyic - dtyiayic; |
| 370 |
dtyibzic = -dtyiczic - dtyiazic; |
| 371 |
dtzibxic = -dtxiczic - dtziaxic; |
| 372 |
dtzibyic = -dtyiczic - dtziayic; |
| 373 |
dtzibzic = -dtziczic - dtziazic; |
| 374 |
dtxibxib = -dtxibxia - dtxibxic; |
| 375 |
dtxibyib = -dtxibyia - dtxibyic; |
| 376 |
dtxibzib = -dtxibzia - dtxibzic; |
| 377 |
dtyibyib = -dtyibyia - dtyibyic; |
| 378 |
dtyibzib = -dtyibzia - dtyibzic; |
| 379 |
dtzibzib = -dtzibzia - dtzibzic; |
| 380 |
|
| 381 |
j = strbnd.isb[i][1]; |
| 382 |
k = strbnd.isb[i][2]; |
| 383 |
dr = 0.0; |
| 384 |
terma = 0.0; |
| 385 |
termc = 0.0; |
| 386 |
if (j > -1) |
| 387 |
{ |
| 388 |
terma = units.stbnunit* strbnd.ksb1[i] ; |
| 389 |
dr = terma * (rab-bonds_ff.bl[j]); |
| 390 |
terma = terma / rab; |
| 391 |
} |
| 392 |
if (k > -1) |
| 393 |
{ |
| 394 |
termc = units.stbnunit* strbnd.ksb2[i]; |
| 395 |
dr = dr + termc*(rcb-bonds_ff.bl[k]); |
| 396 |
termc = termc / rcb; |
| 397 |
} |
| 398 |
ddrdxia = terma * xab; |
| 399 |
ddrdyia = terma * yab; |
| 400 |
ddrdzia = terma * zab; |
| 401 |
ddrdxic = termc * xcb; |
| 402 |
ddrdyic = termc * ycb; |
| 403 |
ddrdzic = termc * zcb; |
| 404 |
ddrdxib = -ddrdxia - ddrdxic; |
| 405 |
ddrdyib = -ddrdyia - ddrdyic; |
| 406 |
ddrdzib = -ddrdzia - ddrdzic; |
| 407 |
|
| 408 |
xab = xab / rab; |
| 409 |
yab = yab / rab; |
| 410 |
zab = zab / rab; |
| 411 |
xcb = xcb / rcb; |
| 412 |
ycb = ycb / rcb; |
| 413 |
zcb = zcb / rcb; |
| 414 |
|
| 415 |
drxiaxia = terma * (1.0-xab*xab); |
| 416 |
drxiayia = -terma * xab*yab; |
| 417 |
drxiazia = -terma * xab*zab; |
| 418 |
dryiayia = terma * (1.0-yab*yab); |
| 419 |
dryiazia = -terma * yab*zab; |
| 420 |
drziazia = terma * (1.0-zab*zab); |
| 421 |
drxicxic = termc * (1.0-xcb*xcb); |
| 422 |
drxicyic = -termc * xcb*ycb; |
| 423 |
drxiczic = -termc * xcb*zcb; |
| 424 |
dryicyic = termc * (1.0-ycb*ycb); |
| 425 |
dryiczic = -termc * ycb*zcb; |
| 426 |
drziczic = termc * (1.0-zcb*zcb); |
| 427 |
drxibxib = drxiaxia + drxicxic; |
| 428 |
drxibyib = drxiayia + drxicyic; |
| 429 |
drxibzib = drxiazia + drxiczic; |
| 430 |
dryibyib = dryiayia + dryicyic; |
| 431 |
dryibzib = dryiazia + dryiczic; |
| 432 |
drzibzib = drziazia + drziczic; |
| 433 |
|
| 434 |
if (ia == iatom) |
| 435 |
{ |
| 436 |
hess.hessx[ia][0] += dt*drxiaxia + dr*dtxiaxia + 2.0*ddtdxia*ddrdxia; |
| 437 |
hess.hessx[ia][1] += dt*drxiayia + dr*dtxiayia + ddtdxia*ddrdyia + ddtdyia*ddrdxia; |
| 438 |
hess.hessx[ia][2] += dt*drxiazia + dr*dtxiazia + ddtdxia*ddrdzia + ddtdzia*ddrdxia; |
| 439 |
hess.hessy[ia][0] += dt*drxiayia + dr*dtxiayia + ddtdyia*ddrdxia + ddtdxia*ddrdyia; |
| 440 |
hess.hessy[ia][1] += dt*dryiayia + dr*dtyiayia + 2.0*ddtdyia*ddrdyia; |
| 441 |
hess.hessy[ia][2] += dt*dryiazia + dr*dtyiazia + ddtdyia*ddrdzia + ddtdzia*ddrdyia; |
| 442 |
hess.hessz[ia][0] += dt*drxiazia + dr*dtxiazia + ddtdzia*ddrdxia + ddtdxia*ddrdzia; |
| 443 |
hess.hessz[ia][1] += dt*dryiazia + dr*dtyiazia + ddtdzia*ddrdyia + ddtdyia*ddrdzia; |
| 444 |
hess.hessz[ia][2] += dt*drziazia + dr*dtziazia + 2.0*ddtdzia*ddrdzia; |
| 445 |
hess.hessx[ib][0] += -dt*drxiaxia + dr*dtxibxia + ddtdxia*ddrdxib + ddtdxib*ddrdxia; |
| 446 |
hess.hessx[ib][1] += -dt*drxiayia + dr*dtxibyia + ddtdxia*ddrdyib + ddtdyib*ddrdxia; |
| 447 |
hess.hessx[ib][2] += -dt*drxiazia + dr*dtxibzia + ddtdxia*ddrdzib + ddtdzib*ddrdxia; |
| 448 |
hess.hessy[ib][0] += -dt*drxiayia + dr*dtyibxia + ddtdyia*ddrdxib + ddtdxib*ddrdyia; |
| 449 |
hess.hessy[ib][1] += -dt*dryiayia + dr*dtyibyia + ddtdyia*ddrdyib + ddtdyib*ddrdyia; |
| 450 |
hess.hessy[ib][2] += -dt*dryiazia + dr*dtyibzia + ddtdyia*ddrdzib + ddtdzib*ddrdyia; |
| 451 |
hess.hessz[ib][0] += -dt*drxiazia + dr*dtzibxia + ddtdzia*ddrdxib + ddtdxib*ddrdzia; |
| 452 |
hess.hessz[ib][1] += -dt*dryiazia + dr*dtzibyia + ddtdzia*ddrdyib + ddtdyib*ddrdzia; |
| 453 |
hess.hessz[ib][2] += -dt*drziazia + dr*dtzibzia + ddtdzia*ddrdzib + ddtdzib*ddrdzia; |
| 454 |
hess.hessx[ic][0] += dr*dtxiaxic + ddtdxia*ddrdxic + ddtdxic*ddrdxia; |
| 455 |
hess.hessx[ic][1] += dr*dtxiayic + ddtdxia*ddrdyic + ddtdyic*ddrdxia; |
| 456 |
hess.hessx[ic][2] += dr*dtxiazic + ddtdxia*ddrdzic + ddtdzic*ddrdxia; |
| 457 |
hess.hessy[ic][0] += dr*dtyiaxic + ddtdyia*ddrdxic + ddtdxic*ddrdyia; |
| 458 |
hess.hessy[ic][1] += dr*dtyiayic + ddtdyia*ddrdyic + ddtdyic*ddrdyia; |
| 459 |
hess.hessy[ic][2] += dr*dtyiazic + ddtdyia*ddrdzic + ddtdzic*ddrdyia; |
| 460 |
hess.hessz[ic][0] += dr*dtziaxic + ddtdzia*ddrdxic + ddtdxic*ddrdzia; |
| 461 |
hess.hessz[ic][1] += dr*dtziayic + ddtdzia*ddrdyic + ddtdyic*ddrdzia; |
| 462 |
hess.hessz[ic][2] += dr*dtziazic + ddtdzia*ddrdzic + ddtdzic*ddrdzia; |
| 463 |
} else if (ib == iatom) |
| 464 |
{ |
| 465 |
hess.hessx[ib][0] += dt*drxibxib + dr*dtxibxib + 2.0*ddtdxib*ddrdxib; |
| 466 |
hess.hessx[ib][1] += dt*drxibyib + dr*dtxibyib + ddtdxib*ddrdyib + ddtdyib*ddrdxib; |
| 467 |
hess.hessx[ib][2] += dt*drxibzib + dr*dtxibzib + ddtdxib*ddrdzib + ddtdzib*ddrdxib; |
| 468 |
hess.hessy[ib][0] += dt*drxibyib + dr*dtxibyib + ddtdyib*ddrdxib + ddtdxib*ddrdyib; |
| 469 |
hess.hessy[ib][1] += dt*dryibyib + dr*dtyibyib + 2.0*ddtdyib*ddrdyib; |
| 470 |
hess.hessy[ib][2] += dt*dryibzib + dr*dtyibzib + ddtdyib*ddrdzib + ddtdzib*ddrdyib; |
| 471 |
hess.hessz[ib][0] += dt*drxibzib + dr*dtxibzib + ddtdzib*ddrdxib + ddtdxib*ddrdzib; |
| 472 |
hess.hessz[ib][1] += dt*dryibzib + dr*dtyibzib + ddtdzib*ddrdyib + ddtdyib*ddrdzib; |
| 473 |
hess.hessz[ib][2] += dt*drzibzib + dr*dtzibzib + 2.0*ddtdzib*ddrdzib; |
| 474 |
hess.hessx[ia][0] += -dt*drxiaxia + dr*dtxibxia + ddtdxib*ddrdxia + ddtdxia*ddrdxib; |
| 475 |
hess.hessx[ia][1] += -dt*drxiayia + dr*dtxibyia + ddtdxib*ddrdyia + ddtdyia*ddrdxib; |
| 476 |
hess.hessx[ia][2] += -dt*drxiazia + dr*dtxibzia + ddtdxib*ddrdzia + ddtdzia*ddrdxib; |
| 477 |
hess.hessy[ia][0] += -dt*drxiayia + dr*dtyibxia + ddtdyib*ddrdxia + ddtdxia*ddrdyib; |
| 478 |
hess.hessy[ia][1] += -dt*dryiayia + dr*dtyibyia + ddtdyib*ddrdyia + ddtdyia*ddrdyib; |
| 479 |
hess.hessy[ia][2] += -dt*dryiazia + dr*dtyibzia + ddtdyib*ddrdzia + ddtdzia*ddrdyib; |
| 480 |
hess.hessz[ia][0] += -dt*drxiazia + dr*dtzibxia + ddtdzib*ddrdxia + ddtdxia*ddrdzib; |
| 481 |
hess.hessz[ia][1] += -dt*dryiazia + dr*dtzibyia + ddtdzib*ddrdyia + ddtdyia*ddrdzib; |
| 482 |
hess.hessz[ia][2] += -dt*drziazia + dr*dtzibzia + ddtdzib*ddrdzia + ddtdzia*ddrdzib; |
| 483 |
hess.hessx[ic][0] += -dt*drxicxic + dr*dtxibxic + ddtdxib*ddrdxic + ddtdxic*ddrdxib; |
| 484 |
hess.hessx[ic][1] += -dt*drxicyic + dr*dtxibyic + ddtdxib*ddrdyic + ddtdyic*ddrdxib; |
| 485 |
hess.hessx[ic][2] += -dt*drxiczic + dr*dtxibzic + ddtdxib*ddrdzic + ddtdzic*ddrdxib; |
| 486 |
hess.hessy[ic][0] += -dt*drxicyic + dr*dtyibxic + ddtdyib*ddrdxic + ddtdxic*ddrdyib; |
| 487 |
hess.hessy[ic][1] += -dt*dryicyic + dr*dtyibyic + ddtdyib*ddrdyic + ddtdyic*ddrdyib; |
| 488 |
hess.hessy[ic][2] += -dt*dryiczic + dr*dtyibzic + ddtdyib*ddrdzic + ddtdzic*ddrdyib; |
| 489 |
hess.hessz[ic][0] += -dt*drxiczic + dr*dtzibxic + ddtdzib*ddrdxic + ddtdxic*ddrdzib; |
| 490 |
hess.hessz[ic][1] += -dt*dryiczic + dr*dtzibyic + ddtdzib*ddrdyic + ddtdyic*ddrdzib; |
| 491 |
hess.hessz[ic][2] += -dt*drziczic + dr*dtzibzic + ddtdzib*ddrdzic + ddtdzic*ddrdzib; |
| 492 |
}else if (ic == iatom) |
| 493 |
{ |
| 494 |
hess.hessx[ic][0] += dt*drxicxic + dr*dtxicxic + 2.0*ddtdxic*ddrdxic; |
| 495 |
hess.hessx[ic][1] += dt*drxicyic + dr*dtxicyic + ddtdxic*ddrdyic + ddtdyic*ddrdxic; |
| 496 |
hess.hessx[ic][2] += dt*drxiczic + dr*dtxiczic + ddtdxic*ddrdzic + ddtdzic*ddrdxic; |
| 497 |
hess.hessy[ic][0] += dt*drxicyic + dr*dtxicyic + ddtdyic*ddrdxic + ddtdxic*ddrdyic; |
| 498 |
hess.hessy[ic][1] += dt*dryicyic + dr*dtyicyic + 2.0*ddtdyic*ddrdyic; |
| 499 |
hess.hessy[ic][2] += dt*dryiczic + dr*dtyiczic + ddtdyic*ddrdzic + ddtdzic*ddrdyic; |
| 500 |
hess.hessz[ic][0] += dt*drxiczic + dr*dtxiczic + ddtdzic*ddrdxic + ddtdxic*ddrdzic; |
| 501 |
hess.hessz[ic][1] += dt*dryiczic + dr*dtyiczic + ddtdzic*ddrdyic + ddtdyic*ddrdzic; |
| 502 |
hess.hessz[ic][2] += dt*drziczic + dr*dtziczic + 2.0*ddtdzic*ddrdzic; |
| 503 |
hess.hessx[ib][0] += -dt*drxicxic + dr*dtxibxic + ddtdxic*ddrdxib + ddtdxib*ddrdxic; |
| 504 |
hess.hessx[ib][1] += -dt*drxicyic + dr*dtxibyic + ddtdxic*ddrdyib + ddtdyib*ddrdxic; |
| 505 |
hess.hessx[ib][2] += -dt*drxiczic + dr*dtxibzic + ddtdxic*ddrdzib + ddtdzib*ddrdxic; |
| 506 |
hess.hessy[ib][0] += -dt*drxicyic + dr*dtyibxic + ddtdyic*ddrdxib + ddtdxib*ddrdyic; |
| 507 |
hess.hessy[ib][1] += -dt*dryicyic + dr*dtyibyic + ddtdyic*ddrdyib + ddtdyib*ddrdyic; |
| 508 |
hess.hessy[ib][2] += -dt*dryiczic + dr*dtyibzic + ddtdyic*ddrdzib + ddtdzib*ddrdyic; |
| 509 |
hess.hessz[ib][0] += -dt*drxiczic + dr*dtzibxic + ddtdzic*ddrdxib + ddtdxib*ddrdzic; |
| 510 |
hess.hessz[ib][1] += -dt*dryiczic + dr*dtzibyic + ddtdzic*ddrdyib + ddtdyib*ddrdzic; |
| 511 |
hess.hessz[ib][2] += -dt*drziczic + dr*dtzibzic + ddtdzic*ddrdzib + ddtdzib*ddrdzic; |
| 512 |
hess.hessx[ia][0] += dr*dtxiaxic + ddtdxic*ddrdxia + ddtdxia*ddrdxic; |
| 513 |
hess.hessx[ia][1] += dr*dtyiaxic + ddtdxic*ddrdyia + ddtdyia*ddrdxic; |
| 514 |
hess.hessx[ia][2] += dr*dtziaxic + ddtdxic*ddrdzia + ddtdzia*ddrdxic; |
| 515 |
hess.hessy[ia][0] += dr*dtxiayic + ddtdyic*ddrdxia + ddtdxia*ddrdyic; |
| 516 |
hess.hessy[ia][1] += dr*dtyiayic + ddtdyic*ddrdyia + ddtdyia*ddrdyic; |
| 517 |
hess.hessy[ia][2] += dr*dtziayic + ddtdyic*ddrdzia + ddtdzia*ddrdyic; |
| 518 |
hess.hessz[ia][0] += dr*dtxiazic + ddtdzic*ddrdxia + ddtdxia*ddrdzic; |
| 519 |
hess.hessz[ia][1] += dr*dtyiazic + ddtdzic*ddrdyia + ddtdyia*ddrdzic; |
| 520 |
hess.hessz[ia][2] += dr*dtziazic + ddtdzic*ddrdzia + ddtdzia*ddrdzic; |
| 521 |
} |
| 522 |
} |
| 523 |
} |
| 524 |
} |
| 525 |
} |
| 526 |
|
