4 * Created by Toshi Nagata on 2005/06/07.
5 * Copyright 2005 Toshi Nagata. All rights reserved.
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation version 2 of the License.
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
18 #include "MDGraphite.h"
24 extern int do_custom_bond_callback(MDArena *arena, Double r, void *procObj, Double *energy, Double *force);
27 s_custom_bond_force(MDArena *arena, Double r, MDCustomBondPar *cp, Double *energy, Double *force)
32 s = exp(-cp->u.morse.a * (r - cp->u.morse.r0));
33 *energy = cp->u.morse.D * (1.0 + s * (s - 2.0));
34 *force = -2.0 * cp->u.morse.a * cp->u.morse.D * s * (1.0 - s);
37 /* retval = do_custom_bond_callback(arena, r, cp->u.proc, energy, force);
39 md_panic(arena, NULL); */
42 *energy = *force = 0.0;
48 s_calc_bond_force(MDArena *arena)
50 Atom *ap = arena->mol->atoms;
51 Int nbonds = arena->mol->nbonds;
52 Int *bonds = arena->mol->bonds;
53 Int *bond_par_i = arena->bond_par_i;
54 Int *custom_bond_par_i = arena->custom_bond_par_i;
55 BondPar *bond_pars = arena->par->bondPars;
56 Double *anbond_r0 = arena->anbond_r0;
57 Double *energies = &arena->energies[kBondIndex];
58 Vector *forces = &arena->forces[kBondIndex * arena->mol->natoms];
60 for (i = 0; i < nbonds; i++, bonds += 2, bond_par_i++) {
62 Double k0, k1, w1, w2, r0;
66 continue; /* Ignore this entry */
67 /* if (bond_uniq != NULL && bond_uniq[i] >= 0)
68 continue; *//* Non-unique bond */
69 if (ap[bonds[0]].mm_exclude || ap[bonds[1]].mm_exclude)
70 continue; /* Skip non-occupied atoms */
71 VecSub(r12, ap[bonds[1]].r, ap[bonds[0]].r);
73 if (custom_bond_par_i != NULL && (idx = custom_bond_par_i[i]) > 0 && idx <= arena->ncustom_bond_pars) {
75 s_custom_bond_force(arena, w1, arena->custom_bond_pars + (idx - 1), &energy, &force);
78 r0 = 0.0; /* To keep complier happy */
80 bp = bond_pars + *bond_par_i;
82 if (anbond_r0 != NULL) {
83 if (anbond_r0[i] > 0.0)
87 k0 = bp->k * w2 * w2; /* Energy */
88 k1 = -2.0 * bp->k * w2 / w1; /* Force / r */
90 VecScaleSelf(r12, k1);
92 VecDec(forces[bonds[0]], r12);
93 VecInc(forces[bonds[1]], r12);
94 if (arena->debug_result && arena->debug_output_level > 1) {
95 fprintf(arena->debug_result, "bond force %d-%d: r=%f, r0=%f, k1=%f, {%f %f %f}\n", bonds[0]+1, bonds[1]+1, w1, r0, k1 * INTERNAL2KCAL, r12.x * INTERNAL2KCAL, r12.y * INTERNAL2KCAL, r12.z * INTERNAL2KCAL);
101 s_calc_angle_force_one(const Vector *r1, const Vector *r2, const Vector *r3, const AnglePar *anp, Double *angle, Double *energy, Double *k1out, Vector *force1, Vector *force2, Vector *force3)
103 Vector r21, r23, f1, f2, v1;
104 Double w1, w2, w3, cost, sint, t, k0, k1;
105 VecSub(r21, *r1, *r2);
106 VecSub(r23, *r3, *r2);
107 w1 = 1.0 / VecLength(r21);
108 w2 = 1.0 / VecLength(r23);
109 VecScaleSelf(r21, w1);
110 VecScaleSelf(r23, w2);
111 cost = VecDot(r21, r23);
112 if (cost > 0.999999 || cost < -0.999999) {
113 /* printf("Cannot handle linear angle %d-%d-%d: skipped.\n", angles[0]+1, angles[1]+1, angles[2]+1); */
114 return -1; /* Cannot handle linear angle */
116 sint = sqrt(1.0 - cost * cost);
117 /* if (sint < 1e-5 && sint > -1e-5)
118 continue; *//* Cannot handle linear angle */
119 t = atan2(sint, cost) - anp->a0;
122 if (sint < 0.1 && sint > -0.1) {
123 /* ---- Method 1 ---- */
124 /* This is slower than method 2, but it is safer when sint is close to 0 */
125 k1 = -2 * anp->k * t;
126 VecCross(v1, r21, r23);
127 VecCross(f1, r21, v1);
128 w3 = w1 * k1 / VecLength(f1);
131 VecScaleSelf(f1, w3);
132 VecCross(f2, v1, r23);
133 w3 = w2 * k1 / VecLength(f2);
136 VecScaleSelf(f2, w3);
138 /* ---- Method 2 ---- */
139 k1 = -2 * anp->k * t / sint;
140 VecScale(f1, r21, cost);
143 VecScaleSelf(f1, w3);
144 VecScale(f2, r23, cost);
147 VecScaleSelf(f2, w3);
153 *angle = t + anp->a0;
156 VecScale(*force2, f1, -1);
161 s_calc_angle_force(MDArena *arena)
163 Atom *ap = arena->mol->atoms;
164 Int nangles = arena->mol->nangles;
165 Int *angles = arena->mol->angles;
166 Int *angle_par_i = arena->angle_par_i;
167 AnglePar *angle_pars = arena->par->anglePars;
168 /* Int *angle_uniq = (arena->nsyms > 0 ? arena->sym_angle_uniq : NULL); */
169 Double *energies = &arena->energies[kAngleIndex];
170 Vector *forces = &arena->forces[kAngleIndex * arena->mol->natoms];
172 for (i = 0; i < nangles; i++, angles += 3, angle_par_i++) {
173 Atom *ap1, *ap2, *ap3;
177 if (*angle_par_i < 0)
178 continue; /* Ignore this entry */
179 ap1 = ATOM_AT_INDEX(ap, angles[0]);
180 ap2 = ATOM_AT_INDEX(ap, angles[1]);
181 ap3 = ATOM_AT_INDEX(ap, angles[2]);
182 /* if (angle_uniq != NULL && angle_uniq[i] >= 0)
183 continue; */ /* Non-unique angle */
184 if (ap1->mm_exclude || ap2->mm_exclude || ap3->mm_exclude)
185 continue; /* Skip non-occupied atoms */
186 if (arena->nalchem_flags > 0) {
187 if (angles[0] < arena->nalchem_flags && angles[1] < arena->nalchem_flags
188 && angles[2] < arena->nalchem_flags
189 && (arena->alchem_flags[angles[0]] | arena->alchem_flags[angles[1]] | arena->alchem_flags[angles[2]]) == 3)
190 continue; /* Interaction between vanishing and appearing groups is ignored */
192 anp = angle_pars + *angle_par_i;
193 if (s_calc_angle_force_one(&(ap1->r), &(ap2->r), &(ap3->r), angle_pars + *angle_par_i, &ang, &en, &k1, &f1, &f2, &f3) != 0)
196 VecInc(forces[angles[0]], f1);
197 VecInc(forces[angles[1]], f2);
198 VecInc(forces[angles[2]], f3);
199 if (arena->debug_result && arena->debug_output_level > 1) {
200 fprintf(arena->debug_result, "angle force %d-%d-%d: a=%f, a0=%f, k1=%f, {%f %f %f}, {%f %f %f}\n", angles[0]+1, angles[1]+1, angles[2]+1, ang*180/PI, anp->a0*180/PI, k1 * INTERNAL2KCAL, -f2.x * INTERNAL2KCAL, -f2.y * INTERNAL2KCAL, -f2.z * INTERNAL2KCAL, f3.x * INTERNAL2KCAL, f3.y * INTERNAL2KCAL, f3.z * INTERNAL2KCAL);
206 s_calc_dihedral_force_one(const Vector *r1, const Vector *r2, const Vector *r3, const Vector *r4, const TorsionPar *tp, Double *phi, Double *energy, Double *k1out, Vector *force1, Vector *force2, Vector *force3, Vector *force4)
208 Vector r21, r32, r43;
210 Double w1, w2, w3, k0, k1;
211 Double cosphi, sinphi;
214 VecSub(r21, *r1, *r2);
215 VecSub(r32, *r2, *r3);
216 VecSub(r43, *r3, *r4);
217 VecCross(v1, r21, r32);
218 VecCross(v2, r32, r43);
219 VecCross(v3, r32, v1);
223 if (w1 < 1e-5 || w2 < 1e-5 || w3 < 1e-5)
224 return -1; /* The dihedral cannot be defined */
228 VecScaleSelf(v1, w1);
229 VecScaleSelf(v2, w2);
230 VecScaleSelf(v3, w3);
232 /* The dihedral angle value */
233 cosphi = VecDot(v1, v2);
234 sinphi = VecDot(v3, v2);
235 *phi = -atan2(sinphi, cosphi);
237 /* Repeat for multiple dihedral terms */
239 for (n = 0; n < tp->mult; n++) {
241 Double phi0 = tp->phi0[n];
242 Int period = tp->period[n];
244 k0 += k * (1 + cos(period * *phi + phi0));
245 k1 -= period * k * sin(period * *phi + phi0);
247 /* A simple quadratic term */
253 k0 += k * phi0 * phi0;
258 if (sinphi < -0.1 || sinphi > 0.1) {
261 VecScale(v4, v1, cosphi);
263 VecScaleSelf(v4, w1);
264 VecScale(v5, v2, cosphi);
266 VecScaleSelf(v5, w2);
268 VecCross(f1, r32, v4);
269 VecCross(f3, v5, r32);
270 VecCross(f2, v4, r21);
271 VecCross(vw0, r43, v5);
273 VecScaleSelf(f1, k1);
274 VecScaleSelf(f2, k1);
275 VecScaleSelf(f3, k1);
278 Vector v6, v7, vw1, vw2;
279 VecScale(v6, v3, sinphi);
281 VecScaleSelf(v6, w3);
282 VecScale(v7, v2, sinphi);
284 VecScaleSelf(v7, w2);
286 VecCross(vw1, v6, r32);
287 VecCross(f1, r32, vw1);
288 VecCross(f3, v7, r32);
289 VecCross(f2, r43, v7);
290 VecCross(vw1, r32, r21);
291 VecCross(vw2, v6, vw1);
293 VecCross(vw1, r32, v6);
294 VecCross(vw2, r21, vw1);
296 VecScaleSelf(f1, k1);
297 VecScaleSelf(f2, k1);
298 VecScaleSelf(f3, k1);
300 /* if (dihedral_uniq != NULL)
301 k0 *= -dihedral_uniq[i]; */
306 VecScale(*force2, f1, -1);
308 VecScale(*force3, f2, -1);
309 VecScale(*force4, f3, -1);
314 s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedrals, Int *dihedral_par_i, TorsionPar *dihedral_pars, Int *dihedral_uniq, Double *energies, Vector *forces)
317 if (arena->mol->nsyms == 0)
318 dihedral_uniq = NULL; /* Ignore the symmetry info */
319 for (i = 0; i < ndihedrals; i++, dihedrals += 4, dihedral_par_i++) {
320 Atom *ap1, *ap2, *ap3, *ap4;
322 Vector f1, f2, f3, f4;
325 if (*dihedral_par_i < 0)
326 continue; /* Ignore this entry */
327 /* if (dihedral_uniq != NULL && dihedral_uniq[i] >= 0)
328 continue; *//* Non-unique dihedral */
329 ap1 = ATOM_AT_INDEX(ap, dihedrals[0]);
330 ap2 = ATOM_AT_INDEX(ap, dihedrals[1]);
331 ap3 = ATOM_AT_INDEX(ap, dihedrals[2]);
332 ap4 = ATOM_AT_INDEX(ap, dihedrals[3]);
333 if (ap1->mm_exclude || ap2->mm_exclude || ap3->mm_exclude || ap4->mm_exclude)
334 continue; /* Skip non-occupied atoms */
335 if (arena->nalchem_flags > 0) {
336 if (dihedrals[0] < arena->nalchem_flags && dihedrals[1] < arena->nalchem_flags
337 && dihedrals[2] < arena->nalchem_flags && dihedrals[3] < arena->nalchem_flags
338 && (arena->alchem_flags[dihedrals[0]] | arena->alchem_flags[dihedrals[1]]
339 | arena->alchem_flags[dihedrals[2]] | arena->alchem_flags[dihedrals[3]]) == 3)
340 continue; /* Interaction between vanishing and appearing groups is ignored */
342 tp = dihedral_pars + *dihedral_par_i;
343 if (s_calc_dihedral_force_one(&(ap1->r), &(ap2->r), &(ap3->r), &(ap4->r), tp, &phi, &en, &k1, &f1, &f2, &f3, &f4) != 0)
347 VecSub(r21, ap[dihedrals[0]].r, ap[dihedrals[1]].r);
348 VecSub(r32, ap[dihedrals[1]].r, ap[dihedrals[2]].r);
349 VecSub(r43, ap[dihedrals[2]].r, ap[dihedrals[3]].r);
350 VecCross(v1, r21, r32);
351 VecCross(v2, r32, r43);
352 VecCross(v3, r32, v1);
356 if (w1 < 1e-5 || w2 < 1e-5 || w3 < 1e-5)
357 continue; // The dihedral cannot be defined
361 VecScaleSelf(v1, w1);
362 VecScaleSelf(v2, w2);
363 VecScaleSelf(v3, w3);
365 // The dihedral angle value
366 cosphi = VecDot(v1, v2);
367 sinphi = VecDot(v3, v2);
368 phi = -atan2(sinphi, cosphi);
370 // Repeat for multiple dihedral terms
372 for (n = 0; n < tp->mult; n++) {
374 Double phi0 = tp->phi0[n];
375 Int period = tp->period[n];
377 k0 += k * (1 + cos(period * phi + phi0));
378 k1 -= period * k * sin(period * phi + phi0);
380 // A simple quadratic term
391 if (sinphi < -0.1 || sinphi > 0.1) {
394 VecScale(v4, v1, cosphi);
396 VecScaleSelf(v4, w1);
397 VecScale(v5, v2, cosphi);
399 VecScaleSelf(v5, w2);
401 VecCross(f1, r32, v4);
402 VecCross(f3, v5, r32);
403 VecCross(f2, v4, r21);
404 VecCross(vw0, r43, v5);
406 VecScaleSelf(f1, k1);
407 VecScaleSelf(f2, k1);
408 VecScaleSelf(f3, k1);
411 Vector v6, v7, vw1, vw2;
412 VecScale(v6, v3, sinphi);
414 VecScaleSelf(v6, w3);
415 VecScale(v7, v2, sinphi);
417 VecScaleSelf(v7, w2);
419 VecCross(vw1, v6, r32);
420 VecCross(f1, r32, vw1);
421 VecCross(f3, v7, r32);
422 VecCross(f2, r43, v7);
423 VecCross(vw1, r32, r21);
424 VecCross(vw2, v6, vw1);
426 VecCross(vw1, r32, v6);
427 VecCross(vw2, r21, vw1);
429 VecScaleSelf(f1, k1);
430 VecScaleSelf(f2, k1);
431 VecScaleSelf(f3, k1);
434 VecInc(forces[dihedrals[0]], f1);
436 VecDec(forces[dihedrals[1]], f1);
438 VecDec(forces[dihedrals[2]], f2);
439 VecDec(forces[dihedrals[3]], f3);
442 VecInc(forces[dihedrals[0]], f1);
443 VecInc(forces[dihedrals[1]], f2);
444 VecInc(forces[dihedrals[2]], f3);
445 VecInc(forces[dihedrals[3]], f4);
446 if (arena->debug_result && arena->debug_output_level > 1) {
447 fprintf(arena->debug_result, "dihedral(improper) force %d-%d-%d-%d: phi=%f, k1=%f, {%f %f %f}, {%f %f %f}, {%f %f %f}\n", dihedrals[0]+1, dihedrals[1]+1, dihedrals[2]+1, dihedrals[3]+1, phi*180/PI, k1 * INTERNAL2KCAL, f1.x * INTERNAL2KCAL, f1.y * INTERNAL2KCAL, f1.z * INTERNAL2KCAL, f2.x * INTERNAL2KCAL, f2.y * INTERNAL2KCAL, f2.z * INTERNAL2KCAL, f3.x * INTERNAL2KCAL, f3.y * INTERNAL2KCAL, f3.z * INTERNAL2KCAL);
453 s_calc_dihedral_force(MDArena *arena)
455 Molecule *mol = arena->mol;
456 Parameter *par = arena->par;
457 Double *energies = &arena->energies[kDihedralIndex];
458 Vector *forces = &arena->forces[kDihedralIndex * arena->mol->natoms];
459 s_calc_dihedral_force_sub(arena, mol->atoms, mol->ndihedrals, mol->dihedrals, arena->dihedral_par_i, par->dihedralPars, NULL/*arena->sym_dihedral_uniq*/, energies, forces);
463 s_calc_improper_force(MDArena *arena)
465 Molecule *mol = arena->mol;
466 Parameter *par = arena->par;
467 Double *energies = &arena->energies[kImproperIndex];
468 Vector *forces = &arena->forces[kImproperIndex * arena->mol->natoms];
469 s_calc_dihedral_force_sub(arena, mol->atoms, mol->nimpropers, mol->impropers, arena->improper_par_i, par->improperPars, NULL/*arena->sym_improper_uniq*/, energies, forces);
472 /* ==================================================================== */
473 /* md_check_verlet_list: only for debugging the verlet list generation */
474 /* ==================================================================== */
476 typedef struct md_check_record {
482 s_md_check_verlet_comparator(const void *ap, const void *bp)
484 Double d = ((const md_check_record *)ap)->len - ((const md_check_record *)bp)->len;
485 return (d < 0 ? -1 : (d > 0 ? 1 : 0));
489 md_check_verlet_list(MDArena *arena)
495 Vector cell_offsets[27];
496 XtalCell *cell = arena->mol->cell;
500 ndx = (arena->periodic_a ? 1 : 0);
501 ndy = (arena->periodic_b ? 1 : 0);
502 ndz = (arena->periodic_c ? 1 : 0);
503 if (cell != NULL && (ndx != 0 || ndy != 0 || ndz != 0)) {
505 for (i = -1; i <= 1; i++) {
506 for (j = -1; j <= 1; j++) {
507 for (k = -1; k <= 1; k++) {
508 VecZero(cell_offsets[nn]);
509 VecScaleInc(cell_offsets[nn], cell->axes[0], i);
510 VecScaleInc(cell_offsets[nn], cell->axes[1], j);
511 VecScaleInc(cell_offsets[nn], cell->axes[2], k);
517 /* Debugging is done with x-y 2-dimensional system */
518 ncr = arena->mol->natoms * (arena->mol->natoms - 1) / 2 * 9;
519 cr = (md_check_record *)malloc(sizeof(md_check_record) * ncr);
521 for (i = 0, api = arena->mol->atoms; i < arena->mol->natoms; i++, api = ATOM_NEXT(api)) {
522 for (j = i + 1, apj = ATOM_AT_INDEX(arena->mol->atoms, j); j < arena->mol->natoms; j++, apj = ATOM_NEXT(apj)) {
524 VecSub(dr, apj->r, api->r);
525 for (dx = -1; dx <= 1; dx++) {
526 for (dy = -1; dy <= 1; dy++) {
528 nn = dx * 9 + dy * 3 + 13;
529 VecInc(dr2, cell_offsets[nn]);
535 cr[k].len = VecLength(dr2);
541 for (k = 0; k < arena->nverlets; k++) {
542 i = arena->verlets[k].n1;
543 j = arena->verlets[k].n2;
546 i = arena->verlets[k].n2;
548 dx = arena->verlets[k].symop.dx;
549 dy = arena->verlets[k].symop.dy;
550 nn = (i * arena->mol->natoms - i * (i + 1) / 2 + (j - i) - 1) * 9 + (dx + 1) * 3 + dy + 1;
551 if (cr[nn].i != i || cr[nn].j != j || cr[nn].dx != dx || cr[nn].dy != dy)
552 fprintf(arena->log_result, "Debug: internal inconsistency\n");
556 mergesort(cr, ncr, sizeof(md_check_record), s_md_check_verlet_comparator);
557 for (i = 0; i < ncr; i++) {
561 len2 = arena->verlets[cr[i].n].length;
562 fprintf(arena->log_result, "Debug: %5d (i=%4d,j=%4d) [dx=%4d, dy=%4d] n=%4d %f %f\n", i, cr[i].i, cr[i].j, cr[i].dx, cr[i].dy, cr[i].n, cr[i].len, len2);
563 if (cr[i].len > arena->pairlist_distance)
568 fprintf(arena->log_result, "Debug: %5d (i=%4d,j=%4d) [dx=%4d, dy=%4d] n=%4d %f %f\n", i, cr[i].i, cr[i].j, cr[i].dx, cr[i].dy, cr[i].n, cr[i].len, arena->verlets[cr[i].n].length);
572 fflush(arena->log_result);
576 /* ==================================================================== */
578 /* Count all lattice points within distance d from point v */
579 /* Ref. "Algorithms for the Shortest and Closest Lattice Vector Problems"
580 Guillaume Hanrot, Xavier Pujol, Damien Stehle
581 Coding and Cryptology; Lecture Notes in Computer Science, 2011 vol. 6639/2011 pp. 159-190
582 DOI: 10.1007/978-3-642-20901-7_10 */
584 s_enum_neighbors(MDArena *arena, Vector v, Double d)
586 Vector bn[3]; /* Base vectors for the cell axes */
587 Double bl[3]; /* Square of the length of base vectors */
588 Double mu[3]; /* Non-diagonal terms for Gram-Schmidt orthogonalization */
589 XtalCell *cell = arena->mol->cell;
590 Int dim; /* Number of periodic axes */
591 Int ijk[3]; /* Renumber the indices. For 1-dimensional case, axes[ijk[0]] is the periodic axis.
592 For 2-dimensional case, axes[ijk[0]] and axes[ijk[1]] are the periodic axes.
593 Only uses exchange of two elements, so that ijk[ijk[i]] == i for i=0,1,2 */
594 Vector axes[3]; /* Renumbered cell axes */
595 Double t[3]; /* Projected vector v on bn[] */
596 Int count; /* Number of results */
597 Int x[3]; /* Lattice point candidate */
602 dim = (arena->periodic_a != 0) + (arena->periodic_b != 0) + (arena->periodic_c != 0);
603 ijk[0] = 0; ijk[1] = 1; ijk[2] = 2;
605 /* Non-periodic case: check whether (0,0,0) is within d or not */
606 if (VecLength(v) < d) {
607 x[0] = x[1] = x[2] = 0;
608 AssignArray(&arena->lattice_offsets, &arena->nlattice_offsets, sizeof(Int) * 3, 0, x);
614 if (arena->periodic_b) {
615 ijk[0] = 1; ijk[1] = 0;
616 } else if (arena->periodic_c) {
617 ijk[0] = 2; ijk[2] = 0;
619 } else if (dim == 2) {
620 if (!arena->periodic_a) {
621 ijk[2] = 0; ijk[0] = 2;
622 } else if (!arena->periodic_b) {
623 ijk[1] = 0; ijk[0] = 1;
626 axes[0] = cell->axes[ijk[0]];
627 axes[1] = cell->axes[ijk[1]];
628 axes[2] = cell->axes[ijk[2]];
630 /* Gram-Schmidt orthogonalization */
632 bl[0] = VecLength2(bn[0]);
633 mu[0] = VecDot(axes[1], bn[0]) / bl[0]; /* bl[0] should be non-zero */
641 VecScaleInc(bn[1], bn[0], -mu[0]);
642 bl[1] = VecLength2(bn[1]);
643 mu[1] = VecDot(axes[2], bn[0]) / bl[0];
644 if (dim == 2 || bl[1] < 1e-10) {
651 mu[2] = VecDot(axes[2], bn[1]) / bl[1];
653 VecScaleInc(bn[2], bn[0], -mu[1]);
654 VecScaleInc(bn[2], bn[1], -mu[2]);
655 bl[2] = VecLength2(bn[2]);
662 /* Project the target vector */
663 t[0] = t[1] = t[2] = 0.0;
664 t[0] = VecDot(v, bn[0]) / bl[0];
666 t[1] = VecDot(v, bn[1]) / bl[1];
668 t[2] = VecDot(v, bn[2]) / bl[2];
674 x[0] = x[1] = x[2] = 0;
675 r[0] = r[1] = r[2] = 0.0;
678 x[i] = ceil(t[i] - d / sqrt(bl[i]));
684 for (j = i + 1; j < dim; j++) {
685 w += x[j] * mu[j + i - 1];
688 w2 = x[i] - t[i] + w;
689 r[i] = w2 * w2 * bl[i];
691 for (j = i; j < dim; j++) {
702 AssignArray(&arena->lattice_offsets, &arena->nlattice_offsets, sizeof(Int) * 3, count, xx);
710 for (j = i + 1; j < dim; j++) {
711 w += x[j] * mu[j + i - 1];
713 x[i] = ceil(t[i] - w - sqrt(w2 / bl[i]));
726 /* Update the Verlet list (for non-bonding interaction) */
728 s_make_verlet_list(MDArena *arena)
730 int i, j, k, n, nn, natoms;
731 int dx, dy, dz, ndx, ndy, ndz;
732 Molecule *mol = arena->mol;
733 Atom *atoms = mol->atoms;
735 MDVerlet *vl = arena->verlets;
736 Vector *vdr = arena->verlets_dr;
737 Double limit, limit2;
739 XtalCell *cell = mol->cell;
740 Parameter *par = arena->par;
743 natoms = mol->natoms;
744 for (i = 0; i < natoms; i++)
747 ndx = (arena->periodic_a ? 1 : 0);
748 ndy = (arena->periodic_b ? 1 : 0);
749 ndz = (arena->periodic_c ? 1 : 0);
750 limit = arena->pairlist_distance;
751 limit2 = limit * limit;
753 use_sym = (arena->mol->nsyms > 0);
754 for (i = 0, api = atoms; i < natoms; i++, api++) {
755 MDExclusion *exinfo = arena->exinfo + i;
756 Int index4 = (exinfo + 1)->index0;
757 Int vdw_idx1 = arena->vdw_par_i[i];
758 Int vdw_idx, vdw_idx2;
759 arena->verlet_i[i] = n;
761 if (api->anchor != NULL)
762 continue; /* Skip pi_anchors */
764 for (j = i, apj = atoms + j; j < natoms; j++, apj++) {
770 /* int dxbase, dybase, dzbase; */
774 if (apj->anchor != NULL)
778 if (api->fix_force < 0 && apj->fix_force < 0)
781 /* Non-occupied atoms */
782 if (api->mm_exclude || apj->mm_exclude)
785 /* If no symmetry expansion, then no interaction with self */
786 if (ndx + ndy + ndz == 0 && i == j)
789 VecSub(rij, apj->r, api->r);
791 /* Calculate the cell offset for the nearest neighbor */
792 if (ndx + ndy + ndz != 0 && apj->periodic_exclude == 0 && cell != NULL) {
793 count = s_enum_neighbors(arena, rij, limit);
795 static Int sZeros[3] = {0, 0, 0};
796 AssignArray(&arena->lattice_offsets, &arena->nlattice_offsets, sizeof(Int) * 3, 0, sZeros);
799 /* Check the specific cutoff table for (i, j) pair */
800 vdw_idx2 = arena->vdw_par_i[j];
801 if (vdw_idx1 < 0 || vdw_idx2 < 0)
802 vdw_idx = par->nvdwPars * par->nvdwPars; /* A null record */
803 else if (vdw_idx1 < vdw_idx2)
804 vdw_idx = vdw_idx1 * par->nvdwPars + vdw_idx2;
806 vdw_idx = vdw_idx2 * par->nvdwPars + vdw_idx1;
807 vdw_cutoff = arena->cutoff;
808 for (k = par->nvdwCutoffPars - 1; k >= 0; k--) {
809 VdwCutoffPar *cr = par->vdwCutoffPars + k;
810 if (((cr->type1 == kAtomTypeWildcard || cr->type1 == api->type) && (cr->type2 == kAtomTypeWildcard || cr->type2 == apj->type)) ||
811 ((cr->type1 == kAtomTypeWildcard || cr->type1 == apj->type) && (cr->type2 == kAtomTypeWildcard || cr->type2 == api->type))) {
812 vdw_cutoff = cr->cutoff;
815 if ((cr->type1 == i && cr->type2 == j) || (cr->type1 == j && cr->type2 == i)) {
816 vdw_cutoff = cr->cutoff;
820 if (vdw_cutoff < 0) {
821 /* A negative value of specific cutoff means "cutoff at r_eq * (-specific_cutoff) */
822 MDVdwCache *cp = &(arena->vdw_cache[vdw_idx]);
823 Double r_eq = pow(cp->par.A / cp->par.B * 2.0, 1.0/6.0);
824 vdw_cutoff = r_eq * (-vdw_cutoff);
827 /* Search for pairs */
828 for (nn = 0; nn < count; nn++) {
830 /* Note the minus signs; s_enum_neighbors finds "the nearest lattice point", so
831 the offsets should be subtracted */
832 dx = -arena->lattice_offsets[nn * 3];
833 dy = -arena->lattice_offsets[nn * 3 + 1];
834 dz = -arena->lattice_offsets[nn * 3 + 2];
835 if (dx == 0 && dy == 0 && dz == 0) {
836 /* Pair within the unit cell */
838 continue; /* No interaction with self */
839 /* Is this pair to be excluded? */
840 for (index0 = exinfo->index0, ip = arena->exlist + index0; index0 < index4; index0++, ip++) {
844 if (index0 < exinfo->index3)
845 continue; /* Special exclusion, 1-2, 1-3 */
847 exflag = 1; /* 1-4 interaction */
848 } else if (apj->periodic_exclude) {
851 VecScaleInc(rij0, cell->axes[0], dx);
852 VecScaleInc(rij0, cell->axes[1], dy);
853 VecScaleInc(rij0, cell->axes[2], dz);
854 /* VecInc(rij0, cell_offsets[nn]); */
858 lenij2 = VecLength2(rij0);
859 if (lenij2 <= limit2) {
861 if (n >= arena->max_nverlets) {
862 arena->max_nverlets += 32;
863 vl = (MDVerlet *)realloc(vl, sizeof(MDVerlet) * arena->max_nverlets);
865 md_panic(arena, "Low memory");
868 vlp->vdw_type = (exflag ? 1 : 0);
874 vlp->symop.alive = (vlp->symop.dx != 0 || vlp->symop.dy != 0 || vlp->symop.dz != 0);
875 vlp->index = vdw_idx;
878 vlp->vdw_cutoff = vdw_cutoff;
879 vlp->length = sqrt(lenij2);
881 } /* end if lenij2 <= limit */
886 arena->verlet_i[natoms] = n;
889 arena->last_verlet_step = arena->step;
891 if (arena->debug_result && arena->debug_output_level > 2) {
892 fprintf(arena->debug_result, "\n Verlet list at step %d\n", arena->step);
893 fprintf(arena->debug_result, " {atom1 atom2 vdw_type (=1 if 1-4 bonded) vdw_index cell_translation}\n");
894 for (i = 0; i < arena->nverlets; i++) {
895 fprintf(arena->debug_result, "{%d %d %d %d %06d}%c", vl[i].n1+1, vl[i].n2+1, vl[i].vdw_type, vl[i].index, (vl[i].symop.dx + 50) * 10000 + (vl[i].symop.dy + 50) * 100 + (vl[i].symop.dz + 50), (i % 4 == 3 ? '\n' : ' '));
898 fprintf(arena->debug_result, "\n");
902 /* Calculate the nonbonded force */
903 /* group_flags[] is used for calculation of pair-specific interaction between
904 two groups of atoms */
906 s_calc_nonbonded_force_sub(MDArena *arena, Double *energies, Double *eenergies, Vector *forces, Vector *eforces, const MDGroupFlags *group_flags_1, const MDGroupFlags *group_flags_2)
910 Double limit, elimit, slimit, vlimit, dielec_r;
911 Double slimit6, slimit12, vlimit6, vlimit12;
912 Double lambda, dlambda;
915 /* MDVdwCache *vdw_cache = arena->vdw_cache; */
916 Atom *atoms = arena->mol->atoms;
917 XtalCell *cell = arena->mol->cell;
918 Vector *eforces_corr;
919 Double *eenergies_corr;
922 if (group_flags_1 == NULL && arena->use_ewald > 0) {
923 eforces_corr = &arena->forces[kESCorrectionIndex * arena->mol->natoms];
924 eenergies_corr = &arena->energies[kESCorrectionIndex];
929 eenergies_corr = NULL;
933 memset(arena->cell_forces, 0, sizeof(Double) * 12);
936 /* Check whether the Verlet list needs update */
937 if (arena->last_verlet_step >= 0) {
938 i = arena->mol->natoms - 1;
939 vdr = arena->verlets_dr + i;
940 if (arena->cutoff > arena->electro_cutoff)
941 limit = (arena->pairlist_distance - arena->cutoff) * 0.5;
943 limit = (arena->pairlist_distance - arena->electro_cutoff) * 0.5;
944 limit = limit * limit;
945 for ( ; i >= 0; i--, vdr--) {
946 if (VecLength2(*vdr) >= limit)
951 s_make_verlet_list(arena);
953 /* Calculate the non-bonded interaction for each pair in the Verlet list */
954 vlimit = arena->cutoff * arena->cutoff;
955 elimit = arena->electro_cutoff * arena->electro_cutoff;
956 slimit = arena->switch_distance * arena->switch_distance;
958 slimit = 0.0; /* No switching */
960 slimit6 = 1/(slimit * slimit * slimit);
961 slimit12 = slimit6 * slimit6;
963 vlimit6 = 1/(vlimit * vlimit * vlimit);
964 vlimit12 = vlimit6 * vlimit6;
965 dielec_r = COULOMBIC / arena->dielectric;
966 /* if (arena->use_ewald)
967 beta2 = arena->ewald_beta * arena->ewald_beta;
969 for (i = arena->nverlets - 1, vl = arena->verlets + i; i >= 0; i--, vl--) {
971 MDVdwCache *vp = arena->vdw_cache + vl->index;
973 Vector rij, fij, rij2, fij2;
974 Double r2, w2, w6, w12;
975 if (group_flags_1 != NULL && group_flags_2 != NULL) {
976 if (!(get_group_flag(group_flags_1, vl->n1) && get_group_flag(group_flags_2, vl->n2))
977 && !(get_group_flag(group_flags_1, vl->n2) && get_group_flag(group_flags_2, vl->n1)))
981 if (arena->nalchem_flags > 0) {
983 if (vl->n1 < arena->nalchem_flags)
984 c1 = arena->alchem_flags[vl->n1];
986 if (vl->n2 < arena->nalchem_flags)
987 c2 = arena->alchem_flags[vl->n2];
989 if ((c1 == 1 && c2 == 2) || (c1 == 2 && c2 == 1))
991 if (c1 == 1 || c2 == 1) {
992 lambda = (1.0 - arena->alchem_lambda);
993 dlambda = -arena->alchem_dlambda;
994 } else if (c1 == 2 || c2 == 2) {
995 lambda = arena->alchem_lambda;
996 dlambda = arena->alchem_dlambda;
1006 if (vl->vdw_type == 1) {
1009 /* vofs = vp->vcut14; */
1013 /* vofs = vp->vcut; */
1015 ap1 = &atoms[vl->n1];
1016 ap2 = &atoms[vl->n2];
1018 if (vl->symop.alive) {
1019 VecScaleInc(rij, cell->axes[0], vl->symop.dx);
1020 VecScaleInc(rij, cell->axes[1], vl->symop.dy);
1021 VecScaleInc(rij, cell->axes[2], vl->symop.dz);
1023 VecDec(rij, ap1->r);
1024 r2 = VecLength2(rij);
1027 limit = vl->vdw_cutoff * vl->vdw_cutoff;
1031 fij.x = fij.y = fij.z = 0.0;
1034 /* Coulombic force */
1035 w12 = ap1->charge * ap2->charge * dielec_r;
1038 if (do_ewald || !arena->use_xplor_shift) {
1039 w2 = 1.0 / sqrt(r2);
1044 k2 = arena->ewald_beta / w2;
1046 k3 = -k1 * (-er + 2 * k2 * exp(-k2 * k2) * PI2R);
1049 /* vofs = -w12 / arena->electro_cutoff;
1050 k0 += vofs; */ /* Discontinue potential at cutoff */
1053 w2 = 1.0 / sqrt(r2);
1054 k0 = r2 / elimit - 1.0;
1057 k1 = (3.0 * r2 / elimit - 2.0) / elimit - 1.0 / r2;
1058 k1 = -w12 * k1 * w2;
1060 if (vl->vdw_type == 1) {
1061 k0 *= arena->scale14_elect;
1062 k1 *= arena->scale14_elect;
1065 VecScale(fij, rij, k3);
1066 *eenergies_corr += k2 / vl->mult;
1067 VecDec(eforces_corr[vl->n1], fij);
1068 VecInc(eforces_corr[vl->n2], fij);
1069 if (arena->debug_result && arena->debug_output_level > 1) {
1070 fprintf(arena->debug_result, "Electrostatic correction force %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", vl->n1+1, vl->n2+1, sqrt(r2), k2/KCAL2INTERNAL, k3*sqrt(r2)/KCAL2INTERNAL, fij.x/KCAL2INTERNAL, fij.y/KCAL2INTERNAL, fij.z/KCAL2INTERNAL);
1073 VecScale(fij, rij, k1);
1074 *eenergies += k0 / vl->mult;
1075 if (eforces != NULL) {
1076 VecDec(eforces[vl->n1], fij);
1077 VecInc(eforces[vl->n2], fij);
1080 if (arena->debug_result && arena->debug_output_level > 1) {
1081 fprintf(arena->debug_result, "nonbonded(electrostatic) force %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", vl->n1+1, vl->n2+1, sqrt(r2), k0/KCAL2INTERNAL, k1*sqrt(r2)/KCAL2INTERNAL, fij.x/KCAL2INTERNAL, fij.y/KCAL2INTERNAL, fij.z/KCAL2INTERNAL);
1085 /* van der Waals force */
1086 /* w2 = 1/(r**2), w6 = 1/(r**6), w12 = 1/(r**12) */
1090 if (slimit == 0 || r2 < slimit) {
1091 k0 = A * w12 - B * w6;
1092 k1 = 6 * w2 * (2.0 * A * w12 - B * w6);
1094 /* Between switch_distance and cutoff: Linear switching function */
1098 k1 = (A * w12 - B * w6) / (arena->cutoff - arena->switch_distance);
1100 k0 = k1 * (arena->cutoff - w2);
1103 if (vl->vdw_type == 1) {
1104 k0 *= arena->scale14_vdw;
1105 k1 *= arena->scale14_vdw;
1109 VecScale(fij, rij, k1);
1110 *energies += k0 * lambda;
1111 if (forces != NULL) {
1112 VecDec(forces[vl->n1], fij);
1113 VecInc(forces[vl->n2], fij);
1117 arena->alchem_energy += k0 * dlambda;
1119 if (arena->debug_result && arena->debug_output_level > 1) {
1120 fprintf(arena->debug_result, "nonbonded(vdw) force %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", vl->n1+1, vl->n2+1, sqrt(r2), k0/KCAL2INTERNAL, k1*sqrt(r2)/KCAL2INTERNAL, fij.x/KCAL2INTERNAL, fij.y/KCAL2INTERNAL, fij.z/KCAL2INTERNAL);
1124 if (arena->is_minimizing && arena->minimize_cell && (ap2->symop.alive || vl->symop.alive)) {
1125 /* Force for changing cell parameters */
1128 Transform da, tr, symtr, symtr2;
1129 if (ap2->symop.alive)
1130 rij2 = ATOM_AT_INDEX(arena->mol->atoms, ap2->symbase)->r;
1132 TransformVec(&rr, cell->rtr, &rij2);
1133 memmove(symtr, SYMMETRY_AT_INDEX(arena->mol->syms, ap2->symop.sym), sizeof(Transform));
1134 symtr[9] += ap2->symop.dx + vl->symop.dx;
1135 symtr[10] += ap2->symop.dy + vl->symop.dy;
1136 symtr[11] += ap2->symop.dz + vl->symop.dz;
1137 TransformMul(symtr2, symtr, cell->rtr);
1138 TransformMul(symtr2, cell->tr, symtr2);
1139 for (j = 0; j < 12; j++) {
1141 if (arena->periodic_a == 0 && (j == 0 || j == 1 || j == 2 || j == 3 || j == 6 || j == 9))
1143 if (arena->periodic_b == 0 && (j == 1 || j == 3 || j == 4 || j == 5 || j == 7 || j == 10))
1145 if (arena->periodic_c == 0 && (j == 2 || j == 5 || j == 6 || j == 7 || j == 8 || j == 11))
1147 memset(da, 0, sizeof(Transform));
1149 TransformMul(tr, symtr2, da);
1150 /* The translation part should be adjusted, because the derivative matrix
1151 has the (3,3) element as 0.0 (instead of 1.0). Thus,
1152 | R v | | R' v' | = | R*R' R*v' | (instead of R*v'+v)
1153 | 0 1 | | 0 0 | | 0 0 |
1154 da * symtr does not have this problem, because the derivative matrix is
1155 multipled from the left. */
1157 tr[10] -= symtr2[10];
1158 tr[11] -= symtr2[11];
1159 TransformMul(da, da, symtr);
1160 for (k = 0; k < 12; k++)
1162 TransformVec(&rrr, da, &rr);
1163 arena->cell_forces[j] += VecDot(fij2, rrr);
1171 /* Calculate correction terms for excluded atom pairs */
1175 Double d, dd, w0, w1, w12;
1176 for (i = 0, api = arena->mol->atoms; i < arena->mol->natoms; i++, api = ATOM_NEXT(api)) {
1177 MDExclusion *exinfo = arena->exinfo + i;
1178 for (k = exinfo->index0; k < exinfo->index3; k++) {
1179 j = arena->exlist[k];
1182 apj = ATOM_AT_INDEX(arena->mol->atoms, j);
1183 VecSub(rij, api->r, apj->r);
1185 w12 = -api->charge * apj->charge * COULOMBIC / arena->dielectric;
1186 dd = arena->ewald_beta * d;
1187 w0 = w12 * erf(dd) / d;
1188 w1 = -(w0 - w12 * 2 * arena->ewald_beta * PI2R * exp(-dd * dd)) / (d * d);
1189 VecScaleSelf(rij, w1);
1190 if (arena->debug_result && arena->debug_output_level > 1) {
1191 fprintf(arena->debug_result, "Electrostatic correction force (excluded) %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", i+1, j+1, d, w0/KCAL2INTERNAL, w0*d/KCAL2INTERNAL, rij.x/KCAL2INTERNAL, rij.y/KCAL2INTERNAL, rij.z/KCAL2INTERNAL);
1193 VecInc(eforces_corr[i], rij);
1194 VecDec(eforces_corr[j], rij);
1195 *eenergies_corr += w0;
1198 if (arena->debug_result && arena->debug_output_level > 0) {
1199 fprintf(arena->debug_result, "Electrostatic correction energy (non-bonded + excluded): %f\n", *eenergies_corr/KCAL2INTERNAL);
1205 s_calc_nonbonded_force(MDArena *arena)
1207 Double *energies = &arena->energies[kVDWIndex];
1208 Double *eenergies = &arena->energies[kElectrostaticIndex];
1209 Vector *forces = &arena->forces[kVDWIndex * arena->mol->natoms];
1210 Vector *eforces = &arena->forces[kElectrostaticIndex * arena->mol->natoms];
1211 s_calc_nonbonded_force_sub(arena, energies, eenergies, forces, eforces, NULL, NULL);
1215 s_calc_auxiliary_force(MDArena *arena)
1217 Double *energies = &arena->energies[kAuxiliaryIndex];
1218 Vector *forces = &arena->forces[kAuxiliaryIndex * arena->mol->natoms];
1219 Atom *atoms = arena->mol->atoms;
1220 int natoms = arena->mol->natoms;
1223 /* Centric force - OBSOLETE */
1224 /* if (arena->centric_potential_force > 0) {
1225 Vector center = arena->mol->atoms[arena->centric_potential_center].r;
1226 Double rin = arena->centric_potential_inner_limit;
1227 Double rout = arena->centric_potential_outer_limit;
1228 Double k = arena->centric_potential_force;
1229 for (i = 0; i < natoms; i++) {
1231 Double w1, w2, k0, k1;
1232 VecSub(r21, atoms[i].r, center);
1233 w2 = VecLength2(r21);
1237 } else if (rout > 0 && w1 > rout) {
1238 k0 = k * (rout - rin) * (w1 + w1 - rout + rin);
1239 k1 = -2.0 * k * (rout - rin) / w1;
1241 k0 = k * (w1 - rin) * (w1 - rin);
1242 k1 = -2.0 * k * (1 - rin / w1);
1245 VecScaleSelf(r21, k1);
1246 VecInc(forces[i], r21);
1248 if (arena->debug_result && arena->debug_output_level > 1) {
1249 fprintf(arena->debug_result, "auxiliary(centric) force %d: r=%f, k1=%f, {%f %f %f}\n", i, w1, k1*w1 * INTERNAL2KCAL, r21.x * INTERNAL2KCAL, r21.y * INTERNAL2KCAL, r21.z * INTERNAL2KCAL);
1256 /* Spherical boundary conditions */
1257 if (arena->spherical_bc_force > 0) {
1258 Vector center = arena->spherical_bc_center;
1259 Double rin = arena->spherical_bc_inner_limit;
1260 Double rout = arena->spherical_bc_outer_limit;
1261 Double k = arena->spherical_bc_force;
1262 for (i = 0; i < natoms; i++) {
1264 Double w1, w2, k0, k1;
1265 if (atoms[i].anchor != NULL)
1267 VecSub(r21, atoms[i].r, center);
1268 w2 = VecLength2(r21);
1272 } else if (rout > 0 && w1 > rout) {
1273 k0 = k * (rout - rin) * (w1 + w1 - rout + rin);
1274 k1 = -2.0 * k * (rout - rin) / w1;
1276 k0 = k * (w1 - rin) * (w1 - rin);
1277 k1 = -2.0 * k * (1 - rin / w1);
1280 VecScaleSelf(r21, k1);
1281 VecInc(forces[i], r21);
1282 if (arena->debug_result && arena->debug_output_level > 1) {
1283 fprintf(arena->debug_result, "auxiliary(spherical BC) force %d: r=%f, k1=%f, {%f %f %f}\n", i, w1, k1*w1 * INTERNAL2KCAL, r21.x * INTERNAL2KCAL, r21.y * INTERNAL2KCAL, r21.z * INTERNAL2KCAL);
1289 if (arena->box_potential_force > 0) {
1290 Double xsize = arena->box_potential_xsize;
1291 Double ysize = arena->box_potential_ysize;
1292 Double zsize = arena->box_potential_zsize;
1293 Double k = arena->box_potential_force;
1294 for (i = 0; i < natoms; i++) {
1295 Vector r = atoms[i].r;
1296 if (atoms[i].anchor != NULL)
1300 else if (r.x < -xsize)
1305 else if (r.y < -ysize)
1310 else if (r.z < -zsize)
1313 *energies += k * (r.x * r.x + r.y * r.y + r.z * r.z);
1314 VecScaleInc(forces[i], r, -2);
1315 if (arena->debug_result && arena->debug_output_level > 1) {
1316 fprintf(arena->debug_result, "auxiliary(box) force %d: {%f %f %f}\n", i, -2*r.x * INTERNAL2KCAL, -2*r.y * INTERNAL2KCAL, -2*r.z * INTERNAL2KCAL);
1322 for (i = j = 0; i < natoms; i++) {
1323 Atom *ap = &atoms[i];
1325 Double w1, w2, k0, k1;
1326 if (ap->fix_force <= 0.0)
1328 VecSub(r21, ap->r, ap->fix_pos);
1329 w2 = VecLength2(r21);
1331 k0 = ap->fix_force * w2;
1332 k1 = -2.0 * ap->fix_force * w1;
1334 VecScaleSelf(r21, k1);
1335 VecInc(forces[i], r21);
1337 if (arena->debug_result && arena->debug_output_level > 1) {
1338 fprintf(arena->debug_result, "auxiliary(fix) force %d: r=%f, k1=%f, {%f %f %f}\n", i, w1, k1*w1 * INTERNAL2KCAL, r21.x * INTERNAL2KCAL, r21.y * INTERNAL2KCAL, r21.z * INTERNAL2KCAL);
1342 /* External forces (used for artificial deformation of the structure) */
1343 if (arena->nexforces > 0) {
1344 for (i = 0; i < arena->nexforces; i++) {
1345 VecInc(forces[i], arena->exforces[i]);
1350 if (arena->graphite != NULL && arena->use_graphite) {
1351 graphite_force(arena->graphite, arena, energies, forces);
1356 s_md_debug_output_positions(MDArena *arena)
1359 if (arena->debug_result == NULL || arena->debug_output_level == 0)
1361 fprintf(arena->debug_result, "\n Atom positions, velocities, and forces at step %d\n", arena->step);
1362 fprintf(arena->debug_result, "%5s %7s %7s %7s %7s %7s %7s %7s %7s %7s\n", "No.", "x", "y", "z", "vx", "vy", "vz", "fx", "fy", "fz");
1363 for (i = 0; i < arena->mol->natoms; i++) {
1364 Atom *ap = &arena->mol->atoms[i];
1365 fprintf(arena->debug_result, "%5d %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", i+1, ap->r.x, ap->r.y, ap->r.z, ap->v.x, ap->v.y, ap->v.z, ap->f.x * INTERNAL2KCAL, ap->f.y * INTERNAL2KCAL, ap->f.z * INTERNAL2KCAL);
1370 calc_force(MDArena *arena)
1374 Molecule *mol = arena->mol;
1375 Int doSurface = 0, doEwald = 0;
1378 natoms = mol->natoms;
1380 /* Clear the energy/force storage */
1381 for (i = 0; i < natoms; i++)
1382 VecZero(mol->atoms[i].f);
1384 memset(arena->energies, 0, sizeof(Double) * kSlowIndex);
1385 memset(arena->forces, 0, sizeof(Vector) * kSlowIndex * natoms);
1386 arena->total_energy = 0.0;
1387 arena->alchem_energy = 0.0;
1389 if (arena->is_minimizing || arena->step == arena->start_step || arena->step % arena->surface_potential_freq == 0) {
1391 arena->energies[kSurfaceIndex] = 0.0;
1392 memset(arena->forces + kSurfaceIndex * natoms, 0, sizeof(Vector) * natoms);
1395 if (arena->is_minimizing || arena->step == arena->start_step || arena->step % arena->ewald_freq == 0) {
1397 arena->energies[kPMEIndex] = 0.0;
1398 memset(arena->forces + kPMEIndex * natoms, 0, sizeof(Vector) * natoms);
1401 s_calc_bond_force(arena);
1402 s_calc_angle_force(arena);
1403 s_calc_dihedral_force(arena);
1404 s_calc_improper_force(arena);
1405 s_calc_nonbonded_force(arena);
1406 s_calc_auxiliary_force(arena);
1408 if (doEwald && arena->use_ewald != 0)
1409 calc_ewald_force(arena);
1411 if (doSurface && arena->probe_radius > 0.0) {
1412 /* if (arena->sphere_points == 0)
1413 calc_surface_force(arena);
1415 calc_surface_force_2(arena); */
1416 calc_surface_force(arena);
1419 /* Distribute forces on pi-anchor atoms to their components */
1420 for (i = 0; i < natoms; i++) {
1423 ap = &(mol->atoms[i]);
1424 if (ap->anchor == NULL)
1426 ip = AtomConnectData(&ap->anchor->connect);
1427 n = ap->anchor->connect.count;
1428 for (k = 0; k < n; k++) {
1429 ff = &arena->forces[ip[k]];
1430 fa = &arena->forces[i];
1431 w = ap->anchor->coeffs[k];
1432 for (j = 0; j < kKineticIndex; j++) {
1433 VecScaleInc(*ff, *fa, w);
1438 fa = &arena->forces[i];
1439 for (j = 0; j < kKineticIndex; j++) {
1445 /* Sum up all partial forces and energies */
1446 arena->total_energy = 0.0;
1447 for (i = 0; i < kKineticIndex; i++)
1448 arena->total_energy += arena->energies[i];
1450 for (i = 0; i < natoms; i++) {
1451 if (mol->atoms[i].mm_exclude)
1453 fa = &mol->atoms[i].f;
1454 ff = &arena->forces[i];
1455 for (j = 0; j < kKineticIndex; j++) {
1461 /* The total (internal) force must be zero */
1466 for (i = 0; i < natoms; i++) {
1467 Vector *fp = &(mol->atoms[i].f);
1470 w = (natoms > 0 ? 1.0 / natoms : 0.0);
1472 for (i = 0; i < natoms; i++) {
1473 Vector *fp = &(mol->atoms[i].f);
1478 s_md_debug_output_positions(arena);
1482 calc_pair_interaction(MDArena *arena, const MDGroupFlags *group_flags_1, const MDGroupFlags *group_flags_2)
1484 Vector *forces, *eforces;
1485 if (arena->pair_forces != NULL) {
1486 forces = arena->pair_forces;
1487 eforces = forces + arena->mol->natoms;
1488 memset(forces, 0, sizeof(Vector) * arena->mol->natoms * 2);
1489 } else forces = eforces = NULL;
1490 arena->pair_energies[0] = arena->pair_energies[1] = 0.0;
1491 s_calc_nonbonded_force_sub(arena, &arena->pair_energies[0], &arena->pair_energies[1], forces, eforces, group_flags_1, group_flags_2);