MolLib/MD/MDCore.c

   1 /*
   2  *  MDCore.c
   3  *
   4  *  Created by Toshi Nagata on 2005/06/06.
   5  *  Copyright 2005 Toshi Nagata. All rights reserved.
   6  *
   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.
  10
  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.
  15  */
  16
  17 #include "MDCore.h"
  18 #include "MDGraphite.h"
  19 #include "MDPressure.h"
  20
  21 #include <stdlib.h>
  22 #include <string.h>
  23 #include <math.h>
  24 #include <time.h>
  25 #include <unistd.h>
  26
  27 #if __WXMSW__
  28 #define ftello(x) ftell(x)
  29 #endif
  30
  31 static char sErrBuf[256];
  32
  33 void
  34 md_panic(MDArena *arena, const char *fmt,...)
  35 {
  36         va_list ap;
  37         jmp_buf *envp;
  38
  39         /*  Clean up the running simulation  */
  40         md_flush_output_files(arena);
  41         arena->is_running = 0;
  42         arena->mol->needsMDRebuild = 1;
  43
  44         va_start(ap, fmt);
  45         vsnprintf(arena->errmsg, sizeof(arena->errmsg), fmt, ap);
  46         va_end(ap);
  47
  48         envp = arena->setjmp_buf;
  49         arena->setjmp_buf = NULL;
  50         if (arena->md_panic_func != NULL)
  51                 (*(arena->md_panic_func))(arena, arena->errmsg);
  52         if (envp != NULL)
  53                 longjmp(*envp, 1);
  54         else {
  55                 fprintf(stderr, "%s\n", sErrBuf);
  56                 exit(1);
  57         }
  58 }
  59
  60 /*  Message output  */
  61 extern int MyAppCallback_showScriptMessage(const char *fmt, ...);
  62
  63 int
  64 md_warning(MDArena *arena, const char *fmt, ...)
  65 {
  66         va_list ap;
  67         char buf[1024];
  68         va_start(ap, fmt);
  69         vsnprintf(buf, sizeof buf, fmt, ap);
  70         va_end(ap);
  71         MyAppCallback_showScriptMessage("%s", buf);
  72         if (arena->log_result != NULL)
  73                 fputs(buf, arena->log_result);
  74         return strlen(buf);
  75 }
  76
  77 int
  78 md_log(MDArena *arena, const char *fmt, ...)
  79 {
  80         va_list ap;
  81         char buf[1024];
  82         if (arena->log_result == NULL)
  83                 return 0;
  84         if (fmt == NULL) {
  85                 fflush(arena->log_result);
  86                 return 0;
  87         }
  88         va_start(ap, fmt);
  89         vsnprintf(buf, sizeof buf, fmt, ap);
  90         va_end(ap);
  91         fputs(buf, arena->log_result);
  92         return strlen(buf);
  93 }
  94
  95 void
  96 md_debug(MDArena *arena, const char *fmt,...)
  97 {
  98         va_list ap;
  99         if (arena->debug_result == NULL || arena->debug_output_level == 0)
 100                 return;
 101         va_start(ap, fmt);
 102         vfprintf(arena->debug_result, fmt, ap);
 103         va_end(ap);
 104         fflush(arena->debug_result);
 105 }
 106
 107 #if DEBUG
 108 #define DEBUG_FIT_COORDINATES 0  /*  set to 1 while debugging md_fit_coordinates */
 109 #endif
 110
 111 /*  Calculate the transform that moves the current coordinates to the reference
 112   coordinates with least displacements. The reference coordinates are given
 113   as one of the snapshots (refno = zero-based)  */
 114 int
 115 md_fit_coordinates(MDArena *arena, Int refno, Double *weights, Transform trans)
 116 {
 117         Atom *ap;
 118         Vector *rvf;
 119         Int natoms, nn;
 120         Vector org1, org2;
 121         Int i, j, k;
 122         Double w, w1;
 123         Mat33 r, q, u;
 124         Double eigen_val[3];
 125         Vector eigen_vec[3];
 126         Vector s[3];
 127         if (arena == NULL || arena->mol == NULL || arena->mol->natoms == 0)
 128                 return 1;  /*  Molecule is empty  */
 129         if (refno < 0 || refno >= arena->nsnapshots)
 130                 return 2;  /*  No such snapshot  */
 131         natoms = arena->mol->natoms;
 132         ap = arena->mol->atoms;
 133         rvf = arena->snapshots[refno]->rvf;
 134
 135         /*  Calculate the weighted center  */
 136         for (j = 0; j < 2; j++) {
 137                 Vector *op = (j == 0 ? &org1 : &org2);
 138                 VecZero(*op);
 139                 w = 0.0;
 140                 for (i = 0; i < natoms; i++) {
 141                         w1 = (weights != NULL ? weights[i] : ap[i].weight);
 142                         if (w1 == 0.0)
 143                                 continue;
 144                         if (j == 0)
 145                                 VecScaleInc(*op, ap[i].r, w1);
 146                         else
 147                                 VecScaleInc(*op, rvf[i * 3], w1);
 148                         w += w1;
 149                 }
 150                 w = 1.0 / w;
 151                 VecScaleSelf(*op, w);
 152         }
 153
 154     /*  R = sum(weight[n] * x[n] * t(y[n]));  */
 155     /*  Matrix to diagonalize = R * tR    */
 156         memset(r, 0, sizeof(Mat33));
 157         memset(q, 0, sizeof(Mat33));
 158         memset(u, 0, sizeof(Mat33));
 159         nn = 0;
 160         for (i = 0; i < natoms; i++) {
 161                 Vector v1, v2;
 162                 w1 = (weights != NULL ? weights[i] : ap[i].weight);
 163                 if (w1 == 0.0)
 164                         continue;
 165                 VecSub(v1, ap[i].r, org1);
 166                 VecSub(v2, rvf[3 * i], org2);
 167                 r[0] += w1 * v2.x * v1.x;
 168                 r[1] += w1 * v2.y * v1.x;
 169                 r[2] += w1 * v2.z * v1.x;
 170                 r[3] += w1 * v2.x * v1.y;
 171                 r[4] += w1 * v2.y * v1.y;
 172                 r[5] += w1 * v2.z * v1.y;
 173                 r[6] += w1 * v2.x * v1.z;
 174                 r[7] += w1 * v2.y * v1.z;
 175                 r[8] += w1 * v2.z * v1.z;
 176 /*              for (j = 0; j < 3; j++) {
 177                         for (k = 0; k < 3; k++) {
 178                                 r[3*j+k] += w1 * VecIndex(&v2, j) * VecIndex(&v1, k);
 179 #if DEBUG_FIT_COORDINATES
 180                                 printf("r(%d,%d) += %.6g * %.6g * %.6g (%.6g)\n", j, k, w1, VecIndex(&v2, j), VecIndex(&v1, k), r[3*j+k]);
 181 #endif
 182                         }
 183                 }
 184 */
 185                 nn++;
 186         }
 187         for (i = 0; i < 9; i++)
 188                 r[i] /= (nn * nn);
 189         for (i = 0; i < 3; i++) {
 190                 for (j = 0; j < 3; j++) {
 191                         for (k = 0; k < 3; k++) {
 192                                 q[j*3+i] += r[k*3+i] * r[k*3+j];
 193                         }
 194                 }
 195         }
 196
 197 #if DEBUG_FIT_COORDINATES
 198         printf("Matrix to diagonalize:\n");
 199         for (i = 0; i < 3; i++) {
 200                 printf("%10.6g %10.6g %10.6g\n", q[i*3], q[i*3+1], q[i*3+2]);
 201         }
 202 #endif
 203
 204         if (MatrixSymDiagonalize(q, eigen_val, eigen_vec) != 0)
 205                 return 3;  /*  Cannot determine the eigenvector  */
 206
 207 #if DEBUG_FIT_COORDINATES
 208         for (i = 0; i < 3; i++) {
 209                 printf("Eigenvalue %d = %.6g\n", i+1, eigen_val[i]);
 210                 printf("Eigenvector %d: %.6g %.6g %.6g\n", i+1, eigen_vec[i].x, eigen_vec[i].y, eigen_vec[i].z);
 211         }
 212 #endif
 213
 214     /*  s[i] = normalize(tR * v[i])  */
 215     /*  U = s0*t(v0) + s1*t(v1) + s2*t(v2)  */
 216         MatrixTranspose(r, r);
 217         for (i = 0; i < 3; i++) {
 218                 MatrixVec(&s[i], r, &eigen_vec[i]);
 219                 w1 = 1.0 / VecLength(s[i]);
 220                 VecScaleSelf(s[i], w1);
 221 #if DEBUG_FIT_COORDINATES
 222                 printf("s[%d] = %.6g %.6g %.6g\n", i+1, s[i].x, s[i].y, s[i].z);
 223 #endif
 224         }
 225 /*      for (i = 0; i < 3; i++) {
 226                 for (j = 0; j < 3; j++) {
 227                         for (k = 0; k < 3; k++) {
 228                                 u[3*i+j] += VecIndex(&s[k], j) * VecIndex(&eigen_vec[k], i);
 229                         }
 230                 }
 231         }
 232 */
 233         for (k = 0; k < 3; k++) {
 234                 u[0] += s[k].x * eigen_vec[k].x;
 235                 u[1] += s[k].x * eigen_vec[k].y;
 236                 u[2] += s[k].x * eigen_vec[k].z;
 237                 u[3] += s[k].y * eigen_vec[k].x;
 238                 u[4] += s[k].y * eigen_vec[k].y;
 239                 u[5] += s[k].y * eigen_vec[k].z;
 240                 u[6] += s[k].z * eigen_vec[k].x;
 241                 u[7] += s[k].z * eigen_vec[k].y;
 242                 u[8] += s[k].z * eigen_vec[k].z;
 243         }
 244
 245         /*  y = U*(x - org1) + org2 = U*x + (org2 - U*org1)  */
 246         MatrixVec(&org1, u, &org1);
 247         VecDec(org2, org1);
 248         for (i = 0; i < 9; i++)
 249                 trans[i] = u[i];
 250         trans[9] = org2.x;
 251         trans[10] = org2.y;
 252         trans[11] = org2.z;
 253
 254         return 0;
 255 }
 256
 257 static void
 258 s_register_missing_parameters(Int **missing, Int *nmissing, Int type, Int t1, Int t2, Int t3, Int t4)
 259 {
 260         Int *mp;
 261         Int i;
 262         for (i = 0, mp = *missing; i < *nmissing; i++, mp += 5) {
 263                 if (mp[0] == type && mp[1] == t1 && mp[2] == t2 && mp[3] == t3 && mp[4] == t4)
 264                         return; /* Already registered */
 265         }
 266         mp = (Int *)AssignArray(missing, nmissing, sizeof(Int)*5, *nmissing, NULL);
 267         if (mp != NULL) {
 268                 mp[0] = type;
 269                 mp[1] = t1;
 270                 mp[2] = t2;
 271                 mp[3] = t3;
 272                 mp[4] = t4;
 273         }
 274 }
 275
 276 /*  Check the bonded atoms and append to results if not already present */
 277 /*  results[] is terminated by -1, hence must be at least (natom+1) size  */
 278 static int
 279 s_check_bonded(Atom *ap, Int *results)
 280 {
 281         Int i, n, *ip;
 282         const Int *cp;
 283         cp = AtomConnects(ap);
 284         for (i = 0; i < ap->nconnects; i++, cp++) {
 285                 n = *cp;
 286                 for (ip = results; *ip >= 0; ip++) {
 287                         if (n == *ip)
 288                                 break;
 289                 }
 290                 if (*ip < 0) {
 291                         *ip++ = n;
 292                         *ip = -1;
 293                 }
 294         }
 295         for (ip = results; *ip >= 0; ip++)
 296                 ;
 297         return ip - results;
 298 }
 299
 300 static void
 301 s_make_exclusion_list(MDArena *arena)
 302 {
 303         Int *results;
 304         Int natoms = arena->mol->natoms;
 305         Atom *atoms = arena->mol->atoms;
 306         MDExclusion *exinfo;
 307         int next_index, i, j;
 308
 309         results = (Int *)calloc(sizeof(Int), natoms + 1);
 310         if (results == NULL)
 311                 md_panic(arena, ERROR_out_of_memory);
 312
 313         if (arena->exlist != NULL) {
 314                 free(arena->exlist);
 315                 arena->exlist = NULL;
 316         }
 317         arena->nexlist = 0;
 318
 319         if (arena->exinfo != NULL)
 320                 free(arena->exinfo);
 321         arena->exinfo = (MDExclusion *)calloc(sizeof(MDExclusion), natoms + 1);
 322         if (arena->exinfo == NULL)
 323                 md_panic(arena, ERROR_out_of_memory);
 324         exinfo = arena->exinfo;
 325
 326         next_index = 0;
 327         for (i = 0; i < natoms; i++) {
 328                 int n;
 329                 exinfo[i].index0 = 0;
 330                 /* special exclusion: only self */
 331                 results[0] = i;
 332                 results[1] = -1;
 333                 exinfo[i].index1 = 1;
 334                 /*  1-2 exclusion (directly bonded)  */
 335                 exinfo[i].index2 = s_check_bonded(&atoms[i], results);
 336                 n = exinfo[i].index2;
 337                 /*  1-3 exclusion: atoms bonded to 1-2 exclusions  */
 338                 for (j = exinfo[i].index1; j < exinfo[i].index2; j++)
 339                         n = s_check_bonded(&atoms[results[j]], results);
 340                 exinfo[i].index3 = n;
 341                 /*  1-4 exclusion: atoms bonded to 1-3 exclusions  */
 342                 for (j = exinfo[i].index2; j < exinfo[i].index3; j++)
 343                         n = s_check_bonded(&atoms[results[j]], results);
 344                 AssignArray(&arena->exlist, &arena->nexlist, sizeof(Int), next_index + n, NULL);
 345                 memcpy(arena->exlist + next_index, results, n * sizeof(Int));
 346                 exinfo[i].index0 += next_index;
 347                 exinfo[i].index1 += next_index;
 348                 exinfo[i].index2 += next_index;
 349                 exinfo[i].index3 += next_index;
 350                 next_index += n;
 351         }
 352         exinfo[natoms].index0 = next_index;  /*  End of exlist  */
 353
 354         free(results);
 355 }
 356
 357 static int
 358 s_lookup_improper_pars(Parameter *par, Int type1, Int type2, Int type3, Int type4)
 359 {
 360         Int idx;
 361         Int t1, t2, t3, t4;
 362         for (idx = par->nimproperPars - 1; idx >= 0; idx--) {
 363                 t1 = par->improperPars[idx].type1;
 364                 t2 = par->improperPars[idx].type2;
 365                 t3 = par->improperPars[idx].type3;
 366                 t4 = par->improperPars[idx].type4;
 367                 if (t1 == -3)
 368                         continue;  /*  Custom parameter  */
 369                 if (t3 >= 0 && t3 != type3)
 370                         continue;
 371                 if ((t1 < 0 || t1 == type1) &&
 372                         (t2 < 0 || t2 == type2) &&
 373                         (t4 < 0 || t4 == type4))
 374                         break;
 375                 if ((t1 < 0 || t1 == type1) &&
 376                         (t2 < 0 || t2 == type4) &&
 377                         (t4 < 0 || t4 == type2))
 378                         break;
 379                 if ((t1 < 0 || t1 == type2) &&
 380                         (t2 < 0 || t2 == type1) &&
 381                         (t4 < 0 || t4 == type4))
 382                         break;
 383                 if ((t1 < 0 || t1 == type2) &&
 384                         (t2 < 0 || t2 == type4) &&
 385                         (t4 < 0 || t4 == type1))
 386                         break;
 387                 if ((t1 < 0 || t1 == type4) &&
 388                         (t2 < 0 || t2 == type1) &&
 389                         (t4 < 0 || t4 == type2))
 390                         break;
 391                 if ((t1 < 0 || t1 == type4) &&
 392                         (t2 < 0 || t2 == type2) &&
 393                         (t4 < 0 || t4 == type1))
 394                         break;
 395         }
 396         return idx;
 397 }
 398
 399 int
 400 md_check_abnormal_bond(MDArena *arena, Molecule *mol, int idx)
 401 {
 402         BondPar *bp;
 403         Atom *ap1, *ap2;
 404         Int idx2;
 405         Double d;
 406         if (mol == NULL)
 407                 mol = arena->mol;
 408         if (arena->par == NULL || idx < 0 || idx >= mol->nbonds || (idx2 = arena->bond_par_i[idx]) < 0 || idx2 >= arena->par->nbondPars)
 409                 return -1;
 410         ap1 = ATOM_AT_INDEX(mol->atoms, mol->bonds[idx * 2]);
 411         ap2 = ATOM_AT_INDEX(mol->atoms, mol->bonds[idx * 2 + 1]);
 412         bp = arena->par->bondPars + idx2;
 413         if (bp->k == 0.0 || bp->r0 == 0.0)
 414                 return 0;
 415         d = MoleculeMeasureBond(mol, &(ap1->r), &(ap2->r));
 416         return (fabs(d / bp->r0 - 1.0) >= 0.2);
 417 }
 418
 419 int
 420 md_check_abnormal_angle(MDArena *arena, Molecule *mol, int idx)
 421 {
 422         AnglePar *anp;
 423         Atom *ap1, *ap2, *ap3;
 424         Int idx2;
 425         Double d;
 426         if (mol == NULL)
 427                 mol = arena->mol;
 428         if (arena->par == NULL || idx < 0 || idx >= mol->nangles || (idx2 = arena->angle_par_i[idx]) < 0 || idx2 >= arena->par->nanglePars)
 429                 return -1;
 430         ap1 = ATOM_AT_INDEX(mol->atoms, mol->angles[idx * 3]);
 431         ap2 = ATOM_AT_INDEX(mol->atoms, mol->angles[idx * 3 + 1]);
 432         ap3 = ATOM_AT_INDEX(mol->atoms, mol->angles[idx * 3 + 2]);
 433         anp = arena->par->anglePars + idx2;
 434         if (anp->k == 0.0 || anp->a0 == 0.0)
 435                 return 0;
 436         d = MoleculeMeasureAngle(mol, &(ap1->r), &(ap2->r), &(ap3->r));
 437         return (fabs(d - anp->a0 * kRad2Deg) >= 20.0);
 438 }
 439
 440 int
 441 md_check_abnormal_dihedral(MDArena *arena, Molecule *mol, int idx)
 442 {
 443         TorsionPar *tp;
 444         Atom *ap1, *ap2, *ap3, *ap4;
 445         Int idx2;
 446         Double d;
 447         if (mol == NULL)
 448                 mol = arena->mol;
 449         if (arena->par == NULL || idx < 0 || idx >= mol->ndihedrals || (idx2 = arena->dihedral_par_i[idx]) < 0 || idx2 >= arena->par->ndihedralPars)
 450                 return -1;
 451         ap1 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[idx * 4]);
 452         ap2 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[idx * 4 + 1]);
 453         ap3 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[idx * 4 + 2]);
 454         ap4 = ATOM_AT_INDEX(mol->atoms, mol->dihedrals[idx * 4 + 3]);
 455         tp = arena->par->dihedralPars + idx2;
 456         if (tp->k[0] == 0.0)
 457                 return 0;
 458         d = MoleculeMeasureDihedral(mol, &(ap1->r), &(ap2->r), &(ap3->r), &(ap4->r));
 459         if (tp->period[0] == 0)
 460                 return (fabs(d - tp->phi0[0] * kRad2Deg) >= 20.0);
 461         else {
 462                 d = (tp->period[0] * (d - tp->phi0[0] * kRad2Deg)) / 360.0 + 0.5;
 463                 d = (d - floor(d) - 0.5) * 360.0; /* map to [-180, 180] */
 464                 return (fabs(d) >= 20.0);
 465         }
 466 }
 467
 468 int
 469 md_check_abnormal_improper(MDArena *arena, Molecule *mol, int idx)
 470 {
 471         TorsionPar *tp;
 472         Atom *ap1, *ap2, *ap3, *ap4;
 473         Int idx2;
 474         Double d;
 475         if (mol == NULL)
 476                 mol = arena->mol;
 477         if (arena->par == NULL || idx < 0 || idx >= mol->ndihedrals || (idx2 = arena->improper_par_i[idx]) < 0 || idx2 >= arena->par->nimproperPars)
 478                 return -1;
 479         ap1 = ATOM_AT_INDEX(mol->atoms, mol->impropers[idx * 4]);
 480         ap2 = ATOM_AT_INDEX(mol->atoms, mol->impropers[idx * 4 + 1]);
 481         ap3 = ATOM_AT_INDEX(mol->atoms, mol->impropers[idx * 4 + 2]);
 482         ap4 = ATOM_AT_INDEX(mol->atoms, mol->impropers[idx * 4 + 3]);
 483         tp = arena->par->improperPars + idx2;
 484         if (tp->k[0] == 0.0)
 485                 return 0;
 486         d = MoleculeMeasureDihedral(mol, &(ap1->r), &(ap2->r), &(ap3->r), &(ap4->r));
 487         if (tp->period[0] == 0)
 488                 return (fabs(d - tp->phi0[0] * kRad2Deg) >= 20.0);
 489         else {
 490                 d = (tp->period[0] * (d - tp->phi0[0] * kRad2Deg)) / 360.0 + 0.5;
 491                 d = (d - floor(d) - 0.5) * 360.0; /* map to [-180, 180] */
 492                 return (fabs(d) >= 20.0);
 493         }
 494 }
 495
 496 static int
 497 s_search_bond(MDArena *arena, int n1, int n2)
 498 {
 499         int i, *ip;
 500         if (n1 < 0 || n1 >= arena->mol->natoms || n2 < 0 || n2 >= arena->mol->natoms)
 501                 return -1;
 502         for (i = 0, ip = arena->mol->bonds; i < arena->mol->nbonds; i++, ip += 2) {
 503                 if ((ip[0] == n1 && ip[1] == n2) || (ip[0] == n2 && ip[1] == n1))
 504                         return i;
 505         }
 506         return -1;
 507 }
 508
 509 static int
 510 s_search_angle(MDArena *arena, int n1, int n2, int n3)
 511 {
 512         int i, *ip;
 513         if (n1 < 0 || n1 >= arena->mol->natoms || n2 < 0 || n2 >= arena->mol->natoms || n3 < 0 || n3 >= arena->mol->natoms)
 514                 return -1;
 515         for (i = 0, ip = arena->mol->angles; i < arena->mol->nangles; i++, ip += 3) {
 516                 if (ip[1] == n2 && ((ip[0] == n1 && ip[2] == n3) || (ip[0] == n3 && ip[2] == n1)))
 517                         return i;
 518         }
 519         return -1;
 520 }
 521
 522 static int
 523 s_search_dihedral(MDArena *arena, int n1, int n2, int n3, int n4)
 524 {
 525         int i, *ip;
 526         if (n1 < 0 || n1 >= arena->mol->natoms || n2 < 0 || n2 >= arena->mol->natoms || n3 < 0 || n3 >= arena->mol->natoms || n4 < 0 || n4 >= arena->mol->natoms)
 527                 return -1;
 528         for (i = 0, ip = arena->mol->dihedrals; i < arena->mol->ndihedrals; i++, ip += 4) {
 529                 if ((ip[0] == n1 && ip[1] == n2 && ip[2] == n3 && ip[3] == n4) || (ip[0] == n4 && ip[1] == n3 && ip[2] == n2 && ip[3] == n1))
 530                         return i;
 531         }
 532         return -1;
 533 }
 534
 535 static int
 536 s_search_improper(MDArena *arena, int n1, int n2, int n3, int n4)
 537 {
 538         int i, *ip;
 539         if (n1 < 0 || n1 >= arena->mol->natoms || n2 < 0 || n2 >= arena->mol->natoms || n3 < 0 || n3 >= arena->mol->natoms || n4 < 0 || n4 >= arena->mol->natoms)
 540                 return -1;
 541         for (i = 0, ip = arena->mol->impropers; i < arena->mol->nimpropers; i++, ip += 4) {
 542                 if ((ip[0] == n1 && ip[1] == n2 && ip[2] == n3 && ip[3] == n4) || (ip[0] == n4 && ip[1] == n3 && ip[2] == n2 && ip[3] == n1))
 543                         return i;
 544         }
 545         return -1;
 546 }
 547
 548 /*  Find vdw parameters and build the in-use list */
 549 static int
 550 s_find_vdw_parameters(MDArena *arena)
 551 {
 552         Int idx, i, j, type1, type2, t1, nmissing;
 553         Double cutoff6, cutoff12;
 554         Parameter *par = arena->par;
 555         Molecule *mol = arena->mol;
 556         Atom *ap;
 557         VdwPar *vp;
 558
 559         if (arena->vdw_par_i != NULL)
 560                 free(arena->vdw_par_i);
 561         arena->vdw_par_i = (Int *)calloc(sizeof(Int), mol->natoms);
 562         if (arena->vdw_par_i == NULL)
 563                 md_panic(arena, ERROR_out_of_memory);
 564
 565         /*  Find the vdw parameter; priority: (1) variant-aware in mol->par, (2) mol->par ignoring variants,
 566             (3) variant->aware in global par, (4) global par ignoring variants  */
 567         for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
 568                 if (mol->par != NULL) {
 569                         vp = ParameterLookupVdwPar(mol->par, ap->type, 0);
 570                 /*      vp = ParameterLookupVdwPar(mol->par, i, kParameterLookupLocal | kParameterLookupNoWildcard);
 571                         if (vp == NULL)
 572                                 vp = ParameterLookupVdwPar(mol->par, ap->type, kParameterLookupLocal | kParameterLookupGlobal); */
 573                         if (vp != NULL) {
 574                                 arena->vdw_par_i[i] = (vp - mol->par->vdwPars) + ATOMS_MAX_NUMBER * 2;
 575                                 continue;
 576                         }
 577                 }
 578                 vp = ParameterLookupVdwPar(gBuiltinParameters, ap->type, 0);
 579                 if (vp != NULL) {
 580                         arena->vdw_par_i[i] = (vp - gBuiltinParameters->vdwPars) + ATOMS_MAX_NUMBER;
 581                         continue;
 582                 }
 583                 /*  Record as missing  */
 584                 vp = ParameterLookupVdwPar(par, ap->type, kParameterLookupMissing | kParameterLookupNoBaseAtomType);
 585                 if (vp == NULL) {
 586                         vp = AssignArray(&par->vdwPars, &par->nvdwPars, sizeof(VdwPar), par->nvdwPars, NULL);
 587                         vp->src = -1;
 588                         vp->type1 = ap->type;
 589                 }
 590                 arena->vdw_par_i[i] = (vp - par->vdwPars);
 591         }
 592         nmissing = par->nvdwPars;
 593
 594         /*  Copy the vdw parameters  */
 595         /*  Atoms with the same vdw parameters point to one common entry  */
 596         /*  Except when the atom appears in one of the pair-specific vdw parameters with
 597             atom index specification; in that case, that atom is given a separate entry  */
 598         for (i = 0, idx = par->nvdwPars; i < mol->natoms; i++) {
 599                 t1 = arena->vdw_par_i[i];
 600                 if (t1 < ATOMS_MAX_NUMBER)
 601                         continue;
 602                 arena->vdw_par_i[i] = idx;
 603                 if (mol->par != NULL) {
 604                         /*  Look up the pair-specific vdw parameters with atom index specification  */
 605                         for (j = mol->par->nvdwpPars - 1; j >= 0; j--) {
 606                                 VdwPairPar *vpp = mol->par->vdwpPars + j;
 607                                 if (vpp->type1 == i || vpp->type2 == i)
 608                                         break;
 609                         }
 610                 } else j = -1;
 611                 if (j < 0) {
 612                         /*  Not found: all other entries with the same vdw parameters share the entry  */
 613                         for (j = i + 1; j < mol->natoms; j++) {
 614                                 if (arena->vdw_par_i[j] == t1)
 615                                         arena->vdw_par_i[j] = idx;
 616                         }
 617                 }
 618                 if (t1 >= ATOMS_MAX_NUMBER * 2)
 619                         vp = mol->par->vdwPars + (t1 - ATOMS_MAX_NUMBER * 2);
 620                 else vp = gBuiltinParameters->vdwPars + (t1 - ATOMS_MAX_NUMBER);
 621                 AssignArray(&par->vdwPars, &par->nvdwPars, sizeof(VdwPar), idx, vp);
 622                 idx++;
 623         }
 624
 625         /*  Debug output  */
 626         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
 627                 fprintf(arena->debug_result, "\n  Atom van der Waals parameters\n");
 628                 fprintf(arena->debug_result, "  No. type vdw_i   sigma    eps    sigma14  eps14\n");
 629                 for (i = 0; i < mol->natoms; i++) {
 630                         char s[8];
 631                         idx = arena->vdw_par_i[i];
 632                         fprintf(arena->debug_result, "%5d %-4s %5d ", i+1, AtomTypeDecodeToString(mol->atoms[i].type, s), idx);
 633                         if (idx >= 0) {
 634                                 Double sigma, eps, sigma2, eps2;
 635                                 sigma = par->vdwPars[idx].A;
 636                                 eps = par->vdwPars[idx].B;
 637                                 if (sigma != 0.0) {
 638                                         eps = 0.25 * eps * eps / sigma;
 639                                         sigma = pow(sigma / eps * 0.25, 1.0/12.0);
 640                                 } else {
 641                                         eps = sigma = 0.0;
 642                                 }
 643                                 sigma2 = par->vdwPars[idx].A14;
 644                                 eps2 = par->vdwPars[idx].B14;
 645                                 if (sigma2 != 0.0) {
 646                                         eps2 = 0.25 * eps2 * eps2 / sigma2;
 647                                         sigma2 = pow(sigma2 / eps2 * 0.25, 1.0/12.0);
 648                                 } else {
 649                                         eps2 = sigma2 = 0.0;
 650                                 }
 651                                 fprintf(arena->debug_result, " %7.3f %7.3f %7.3f %7.3f\n", sigma, eps*INTERNAL2KCAL, sigma2, eps2*INTERNAL2KCAL);
 652                         } else {
 653                                 fprintf(arena->debug_result, " (not available)\n");
 654                         }
 655                 }
 656         }
 657
 658         /*  Build cache  */
 659         if (arena->vdw_cache != NULL)
 660                 free(arena->vdw_cache);
 661         arena->vdw_cache = (MDVdwCache *)calloc(sizeof(MDVdwCache), par->nvdwPars * par->nvdwPars + 1);
 662         if (arena->vdw_cache == NULL)
 663                 md_panic(arena, ERROR_out_of_memory);
 664         cutoff6 = arena->cutoff * arena->cutoff;
 665         cutoff6 = 1.0 / (cutoff6 * cutoff6 * cutoff6);
 666         cutoff12 = cutoff6 * cutoff6;
 667         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
 668                 fprintf(arena->debug_result, "\n  van der Waals parameter cache table\n");
 669                 fprintf(arena->debug_result, " idx1  idx2  type1 type2  sigma    eps    sigma14  eps14\n");
 670         }
 671         for (i = 0; i < par->nvdwPars; i++) {
 672         /*      t1 = arena->vdw_par_i[i];
 673                 type1 = par->vdwPars[t1].type1; */
 674                 type1 = par->vdwPars[i].type1;
 675                 for (j = i; j < par->nvdwPars; j++) {
 676                         VdwPairPar newpar, *vpp;
 677                         Double sigma1, sigma2, eps1, eps2;
 678                         vpp = NULL;
 679                 /*      t2 = arena->vdw_par_i[j];
 680                         type2 = par->vdwPars[t2].type1; */
 681                         type2 = par->vdwPars[j].type1;  /*  Not type2  */
 682                         /*  Look up the pair-specific van der Waals parameters  */
 683                         if (mol->par != NULL) {
 684                                 vpp = ParameterLookupVdwPairPar(mol->par, type1, type2, 0);
 685                         }
 686                         if (vpp == NULL)
 687                                 vpp = ParameterLookupVdwPairPar(gBuiltinParameters, type1, type2, 0);
 688                         if (vpp != NULL) {
 689                                 newpar = *vpp;
 690                         } else {
 691                                 /*  Create a pair-specific VdwPar for this pair  */
 692                                 VdwPar *p1, *p2;
 693                                 p1 = par->vdwPars + i;
 694                                 p2 = par->vdwPars + j;
 695                                 if (p1->A < 1e-15 && p1->A > -1e-15) {
 696                                         sigma1 = 1.0;
 697                                         eps1 = 0.0;
 698                                 } else {
 699                                         eps1 = 0.25 * p1->B * p1->B / p1->A;
 700                                         sigma1 = pow(p1->B * 0.25 / eps1, 1.0/6.0);
 701                                 }
 702                                 if (p2->A < 1e-15 && p2->A > -1e-15) {
 703                                         sigma2 = 1.0;
 704                                         eps2 = 0.0;
 705                                 } else {
 706                                         eps2 = 0.25 * p2->B * p2->B / p2->A;
 707                                         sigma2 = pow(p2->B * 0.25 / eps2, 1.0/6.0);
 708                                 }
 709                                 sigma1 = (sigma1 + sigma2) * 0.5;
 710                                 eps1 = sqrt(eps1 * eps2);
 711                                 sigma1 = sigma1 * sigma1 * sigma1;
 712                                 sigma1 = sigma1 * sigma1;
 713                                 newpar.B = 4.0 * sigma1 * eps1;
 714                                 newpar.A = newpar.B * sigma1;
 715                                 if (p1->A14 < 1e-15 && p1->A14 > -1e-15) {
 716                                         sigma1 = 1.0;
 717                                         eps1 = 0.0;
 718                                 } else {
 719                                         eps1 = 0.25 * p1->B14 * p1->B14 / p1->A14;
 720                                         sigma1 = pow(p1->B14 * 0.25 / eps1, 1.0/6.0);
 721                                 }
 722                                 if (p2->A14 < 1e-15 && p2->A14 > -1e-15) {
 723                                         sigma2 = 1.0;
 724                                         eps2 = 0.0;
 725                                 } else {
 726                                         eps2 = 0.25 * p2->B14 * p2->B14 / p2->A14;
 727                                         sigma2 = pow(p2->B14 * 0.25 / eps2, 1.0/6.0);
 728                                 }
 729                                 sigma1 = (sigma1 + sigma2) * 0.5;
 730                                 eps1 = sqrt(eps1 * eps2);
 731                                 sigma1 = sigma1 * sigma1 * sigma1;
 732                                 sigma1 = sigma1 * sigma1;
 733                                 newpar.B14 = 4.0 * sigma1 * eps1;
 734                                 newpar.A14 = newpar.B14 * sigma1;
 735                                 newpar.type1 = type1;
 736                                 newpar.type2 = type2;
 737                         }
 738                         idx = i * par->nvdwPars + j;
 739                         arena->vdw_cache[idx].par = newpar;
 740                         /*  Cache the value at the cutoff distance  */
 741                         arena->vdw_cache[idx].vcut = newpar.A * cutoff12 - newpar.B * cutoff6;
 742                         arena->vdw_cache[idx].vcut14 = newpar.A14 * cutoff12 - newpar.B * cutoff6;
 743                         arena->vdw_cache[j * par->nvdwPars + i] = arena->vdw_cache[idx];
 744                         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
 745                                 char s1[8], s2[8];
 746                                 eps1 = 0.25 * newpar.B * newpar.B / newpar.A;
 747                                 sigma1 = pow(newpar.B * 0.25 / eps1, 1.0/6.0);
 748                                 eps2 = 0.25 * newpar.B14 * newpar.B14 / newpar.A14;
 749                                 sigma2 = pow(newpar.B14 * 0.25 / eps2, 1.0/6.0);
 750                                 fprintf(arena->debug_result, "%5d %5d %-4s  %-4s  %7.3f %7.3f %7.3f %7.3f\n", idx, i, AtomTypeDecodeToString(type1, s1), AtomTypeDecodeToString(type2, s2), sigma1, eps1*INTERNAL2KCAL, sigma2, eps2*INTERNAL2KCAL);
 751                         }
 752                 } /* end loop j */
 753         } /* end loop i */
 754         arena->nmissing += nmissing;
 755         return nmissing;
 756 }
 757
 758 /*  Find bond parameters  */
 759 static int
 760 s_find_bond_parameters(MDArena *arena)
 761 {
 762         Int i, j, idx, type1, type2, t1, nmissing = 0;
 763         Molecule *mol = arena->mol;
 764         Parameter *par = arena->par;
 765         BondPar *bp;
 766
 767         if (mol->nbonds > 0) {
 768                 if (arena->bond_par_i != NULL)
 769                         free(arena->bond_par_i);
 770                 arena->bond_par_i = (Int *)calloc(sizeof(Int), mol->nbonds);
 771                 if (arena->bond_par_i == NULL)
 772                         md_panic(arena, ERROR_out_of_memory);
 773                 if (arena->anbond_r0 != NULL)
 774                         free(arena->anbond_r0);
 775                 arena->anbond_r0 = (Double *)calloc(sizeof(Double), mol->nbonds);
 776                 if (arena->anbond_r0 == NULL)
 777                         md_panic(arena, ERROR_out_of_memory);
 778
 779                 /*  Find the bond parameter; priority: (1) atom index specific in mol->par, (2)
 780                  atom type specific in mol->par, (3) atom type specific in built-in  */
 781                 for (i = 0; i < mol->nbonds; i++) {
 782                         Int i1, i2;
 783                         Atom *ap1, *ap2;
 784                         i1 = mol->bonds[i * 2];
 785                         i2 = mol->bonds[i * 2 + 1];
 786                         ap1 = ATOM_AT_INDEX(mol->atoms, i1);
 787                         ap2 = ATOM_AT_INDEX(mol->atoms, i2);
 788                         type1 = ap1->type;
 789                         type2 = ap2->type;
 790                         bp = NULL;
 791                         if (mol->par != NULL) {
 792                                 bp = ParameterLookupBondPar(mol->par, type1, type2, i1, i2, 0);
 793                                 if (bp != NULL) {
 794                                         arena->bond_par_i[i] = (bp - mol->par->bondPars) + ATOMS_MAX_NUMBER * 2;
 795                                         continue;
 796                                 }
 797                         }
 798                         bp = ParameterLookupBondPar(gBuiltinParameters, type1, type2, -1, -1, 0);
 799                         if (bp != NULL) {
 800                                 arena->bond_par_i[i] = (bp - gBuiltinParameters->bondPars) + ATOMS_MAX_NUMBER;
 801                                 continue;
 802                         }
 803                         bp = ParameterLookupBondPar(par, type1, type2, -1, -1, kParameterLookupMissing | kParameterLookupNoBaseAtomType);
 804                         if (bp == NULL) {
 805                                 /*  Record as missing  */
 806                                 bp = AssignArray(&par->bondPars, &par->nbondPars, sizeof(BondPar), par->nbondPars, NULL);
 807                                 bp->src = -1;
 808                                 bp->type1 = type1;
 809                                 bp->type2 = type2;
 810                         }
 811                         arena->bond_par_i[i] = (bp - par->bondPars);
 812                 }
 813                 nmissing = par->nbondPars;
 814
 815                 /*  Copy the bond parameters  */
 816                 for (i = 0, idx = par->nbondPars; i < mol->nbonds; i++) {
 817                         t1 = arena->bond_par_i[i];
 818                         if (t1 < ATOMS_MAX_NUMBER)
 819                                 continue;
 820                         arena->bond_par_i[i] = idx;
 821                         for (j = i + 1; j < mol->nbonds; j++) {
 822                                 if (arena->bond_par_i[j] == t1)
 823                                         arena->bond_par_i[j] = idx;
 824                         }
 825                         if (t1 >= ATOMS_MAX_NUMBER * 2)
 826                                 bp = mol->par->bondPars + (t1 - ATOMS_MAX_NUMBER * 2);
 827                         else bp = gBuiltinParameters->bondPars + (t1 - ATOMS_MAX_NUMBER);
 828                         AssignArray(&par->bondPars, &par->nbondPars, sizeof(BondPar), idx, bp);
 829                         idx++;
 830                 }
 831         }
 832         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
 833                 char s1[8], s2[8];
 834                 fprintf(arena->debug_result, "\n  Bond parameters\n");
 835                 fprintf(arena->debug_result, "  No. atom1 atom2 type1 type2     r0      k\n");
 836                 for (i = 0; i < mol->nbonds; i++) {
 837                         idx = arena->bond_par_i[i];
 838                         fprintf(arena->debug_result, "%5d %5d %5d ", i+1, mol->bonds[i*2]+1, mol->bonds[i*2+1]+1);
 839                         if (idx < 0) {
 840                                 fprintf(arena->debug_result, "<< missing >>\n");
 841                                 continue;
 842                         }
 843                         fprintf(arena->debug_result, " %-4s  %-4s %7.3f %7.3f\n", AtomTypeDecodeToString(par->bondPars[idx].type1, s1), AtomTypeDecodeToString(par->bondPars[idx].type2, s2), par->bondPars[idx].r0, par->bondPars[idx].k*INTERNAL2KCAL);
 844                 }
 845         }
 846         arena->nmissing += nmissing;
 847         return nmissing;
 848 }
 849
 850 /*  Find angle parameters  */
 851 static int
 852 s_find_angle_parameters(MDArena *arena)
 853 {
 854         Int i, j, idx, type1, type2, type3, t1, nmissing = 0;
 855         Molecule *mol = arena->mol;
 856         Parameter *par = arena->par;
 857         AnglePar *ap;
 858
 859         if (mol->nangles > 0) {
 860                 if (arena->angle_par_i != NULL)
 861                         free(arena->angle_par_i);
 862                 arena->angle_par_i = (Int *)calloc(sizeof(Int), mol->nangles);
 863                 if (arena->angle_par_i == NULL)
 864                         md_panic(arena, ERROR_out_of_memory);
 865
 866                 /*  Find the angle parameter; priority: (1) atom index specific in mol->par, (2)
 867                  atom type specific in mol->par, (3) atom type specific in built-in  */
 868                 for (i = 0; i < mol->nangles; i++) {
 869                         Int i1, i2, i3;
 870                         i1 = mol->angles[i * 3];
 871                         i2 = mol->angles[i * 3 + 1];
 872                         i3 = mol->angles[i * 3 + 2];
 873                         type1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
 874                         type2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
 875                         type3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
 876                         if (mol->par != NULL) {
 877                                 ap = ParameterLookupAnglePar(mol->par, type1, type2, type3, i1, i2, i3, 0);
 878                                 if (ap != NULL) {
 879                                         arena->angle_par_i[i] = (ap - mol->par->anglePars) + ATOMS_MAX_NUMBER * 2;
 880                                         continue;
 881                                 }
 882                         }
 883                         ap = ParameterLookupAnglePar(gBuiltinParameters, type1, type2, type3, -1, -1, -1, 0);
 884                         if (ap != NULL) {
 885                                 arena->angle_par_i[i] = (ap - gBuiltinParameters->anglePars) + ATOMS_MAX_NUMBER;
 886                                 continue;
 887                         }
 888                         /*  Record as missing  */
 889                         ap = ParameterLookupAnglePar(par, type1, type2, type3, -1, -1, -1, kParameterLookupMissing | kParameterLookupNoBaseAtomType);
 890                         if (ap == NULL) {
 891                                 ap = AssignArray(&par->anglePars, &par->nanglePars, sizeof(AnglePar), par->nanglePars, NULL);
 892                                 ap->src = -1;
 893                                 ap->type1 = type1;
 894                                 ap->type2 = type2;
 895                                 ap->type3 = type3;
 896                         }
 897                         arena->angle_par_i[i] = (ap - par->anglePars);
 898                 }
 899                 nmissing = par->nanglePars;
 900
 901                 /*  Copy the angle parameters  */
 902                 for (i = 0, idx = par->nanglePars; i < mol->nangles; i++) {
 903                         t1 = arena->angle_par_i[i];
 904                         if (t1 < ATOMS_MAX_NUMBER)
 905                                 continue;
 906                         arena->angle_par_i[i] = idx;
 907                         for (j = i + 1; j < mol->nangles; j++) {
 908                                 if (arena->angle_par_i[j] == t1)
 909                                         arena->angle_par_i[j] = idx;
 910                         }
 911                         if (t1 >= ATOMS_MAX_NUMBER * 2)
 912                                 ap = mol->par->anglePars + (t1 - ATOMS_MAX_NUMBER * 2);
 913                         else ap = gBuiltinParameters->anglePars + (t1 - ATOMS_MAX_NUMBER);
 914                         AssignArray(&par->anglePars, &par->nanglePars, sizeof(AnglePar), idx, ap);
 915                         idx++;
 916                 }
 917         }
 918
 919         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
 920                 char s1[8], s2[8], s3[8];
 921                 fprintf(arena->debug_result, "\n  Angle parameters\n");
 922                 fprintf(arena->debug_result, "  No. atom1 atom2 atom3 type1 type2 type3      a0      k\n");
 923                 for (i = 0; i < mol->nangles; i++) {
 924                         idx = arena->angle_par_i[i];
 925                         fprintf(arena->debug_result, "%5d %5d %5d %5d ", i+1, mol->angles[i*3]+1, mol->angles[i*3+1]+1, mol->angles[i*3+2]+1);
 926                         if (idx < 0) {
 927                                 fprintf(arena->debug_result, "<< missing >>\n");
 928                                 continue;
 929                         }
 930                         fprintf(arena->debug_result, " %-4s  %-4s  %-4s %7.3f %7.3f\n",  AtomTypeDecodeToString(par->anglePars[idx].type1, s1), AtomTypeDecodeToString(par->anglePars[idx].type2, s2), AtomTypeDecodeToString(par->anglePars[idx].type3, s3), par->anglePars[idx].a0 * 180.0 / 3.1415927, par->anglePars[idx].k*INTERNAL2KCAL);
 931                 }
 932         }
 933         arena->nmissing += nmissing;
 934         return nmissing;
 935 }
 936
 937 /*  Find dihedral parameters  */
 938 static int
 939 s_find_dihedral_parameters(MDArena *arena)
 940 {
 941         Int i, j, idx, type1, type2, type3, type4, t1, nmissing = 0;
 942         Molecule *mol = arena->mol;
 943         Parameter *par = arena->par;
 944         TorsionPar *tp;
 945
 946         if (mol->ndihedrals > 0) {
 947                 if (arena->dihedral_par_i != NULL)
 948                         free(arena->dihedral_par_i);
 949                 arena->dihedral_par_i = (Int *)calloc(sizeof(Int), mol->ndihedrals);
 950                 if (arena->dihedral_par_i == NULL)
 951                         md_panic(arena, ERROR_out_of_memory);
 952
 953                 /*  Find the dihedral parameter; priority: (1) atom index specific in mol->par, (2)
 954                  atom type specific in mol->par, (3) atom type specific in built-in  */
 955                 for (i = 0; i < mol->ndihedrals; i++) {
 956                         Int i1, i2, i3, i4;
 957                         i1 = mol->dihedrals[i * 4];
 958                         i2 = mol->dihedrals[i * 4 + 1];
 959                         i3 = mol->dihedrals[i * 4 + 2];
 960                         i4 = mol->dihedrals[i * 4 + 3];
 961                         type1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
 962                         type2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
 963                         type3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
 964                         type4 = ATOM_AT_INDEX(mol->atoms, i4)->type;
 965                         if (mol->par != NULL) {
 966                                 tp = ParameterLookupDihedralPar(mol->par, type1, type2, type3, type4, i1, i2, i3, i4, 0);
 967                                 if (tp != NULL) {
 968                                         arena->dihedral_par_i[i] = (tp - mol->par->dihedralPars) + ATOMS_MAX_NUMBER * 2;
 969                                         continue;
 970                                 }
 971                         }
 972                         tp = ParameterLookupDihedralPar(gBuiltinParameters, type1, type2, type3, type4, -1, -1, -1, -1, 0);
 973                         if (tp != NULL) {
 974                                 arena->dihedral_par_i[i] = (tp - gBuiltinParameters->dihedralPars) + ATOMS_MAX_NUMBER;
 975                                 continue;
 976                         }
 977                         /*  Record as missing  */
 978                         tp = ParameterLookupDihedralPar(par, type1, type2, type3, type4, -1, -1, -1, -1, kParameterLookupMissing | kParameterLookupNoBaseAtomType);
 979                         if (tp == NULL) {
 980                                 tp = AssignArray(&par->dihedralPars, &par->ndihedralPars, sizeof(TorsionPar), par->ndihedralPars, NULL);
 981                                 tp->src = -1;
 982                                 tp->type1 = type1;
 983                                 tp->type2 = type2;
 984                                 tp->type3 = type3;
 985                                 tp->type4 = type4;
 986                         }
 987                         arena->dihedral_par_i[i] = (tp - par->dihedralPars);
 988                 }
 989                 nmissing = par->ndihedralPars;
 990
 991                 /*  Copy the dihedral parameters  */
 992                 for (i = 0, idx = par->ndihedralPars; i < mol->ndihedrals; i++) {
 993                         t1 = arena->dihedral_par_i[i];
 994                         if (t1 < ATOMS_MAX_NUMBER)
 995                                 continue;
 996                         arena->dihedral_par_i[i] = idx;
 997                         for (j = i + 1; j < mol->ndihedrals; j++) {
 998                                 if (arena->dihedral_par_i[j] == t1)
 999                                         arena->dihedral_par_i[j] = idx;
1000                         }
1001                         if (t1 >= ATOMS_MAX_NUMBER * 2)
1002                                 tp = mol->par->dihedralPars + (t1 - ATOMS_MAX_NUMBER * 2);
1003                         else tp = gBuiltinParameters->dihedralPars + (t1 - ATOMS_MAX_NUMBER);
1004                         AssignArray(&par->dihedralPars, &par->ndihedralPars, sizeof(TorsionPar), idx, tp);
1005                         idx++;
1006                 }
1007         }
1008
1009         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
1010                 Int j;
1011                 char s1[8], s2[8], s3[8], s4[8];
1012                 fprintf(arena->debug_result, "\n  Dihedral parameters\n");
1013                 fprintf(arena->debug_result, "  No. atom1 atom2 atom3 atom4 type1 type2 type3 type4 mult  phi0      k     per\n");
1014                 for (i = 0; i < mol->ndihedrals; i++) {
1015                         idx = arena->dihedral_par_i[i];
1016                         fprintf(arena->debug_result, "%5d %5d %5d %5d %5d ", i+1, mol->dihedrals[i*4]+1, mol->dihedrals[i*4+1]+1, mol->dihedrals[i*4+2]+1, mol->dihedrals[i*4+3]+1);
1017                         if (idx < 0) {
1018                                 fprintf(arena->debug_result, "<< missing >>\n");
1019                                 continue;
1020                         }
1021                         fprintf(arena->debug_result, "%-4s  %-4s  %-4s  %-4s  %2d ", AtomTypeDecodeToString(par->dihedralPars[idx].type1, s1), AtomTypeDecodeToString(par->dihedralPars[idx].type2, s2), AtomTypeDecodeToString(par->dihedralPars[idx].type3, s3), AtomTypeDecodeToString(par->dihedralPars[idx].type4, s4), par->dihedralPars[idx].mult);
1022                         for (j = 0; j < par->dihedralPars[idx].mult; j++) {
1023                                 fprintf(arena->debug_result, "%7.3f %7.3f %1d ", par->dihedralPars[idx].phi0[j]*180/PI, par->dihedralPars[idx].k[j]*INTERNAL2KCAL, par->dihedralPars[idx].period[j]);
1024                         }
1025                         fprintf(arena->debug_result, "\n");
1026                 }
1027         }
1028         arena->nmissing += nmissing;
1029         return nmissing;
1030 }
1031
1032 /*  Find improper parameters  */
1033 static int
1034 s_find_improper_parameters(MDArena *arena)
1035 {
1036         Int i, j, idx, type1, type2, type3, type4, t1, nmissing = 0;
1037         Molecule *mol = arena->mol;
1038         Parameter *par = arena->par;
1039         TorsionPar *tp;
1040
1041         if (mol->nimpropers > 0) {
1042                 if (arena->improper_par_i != NULL)
1043                         free(arena->improper_par_i);
1044                 arena->improper_par_i = (Int *)calloc(sizeof(Int), mol->nimpropers);
1045                 if (arena->improper_par_i == NULL)
1046                         md_panic(arena, ERROR_out_of_memory);
1047
1048                 /*  Find the improper parameter; priority: (1) atom index specific in mol->par, (2)
1049                  atom type specific in mol->par, (3) atom type specific in built-in  */
1050                 for (i = 0; i < mol->nimpropers; i++) {
1051                         Int i1, i2, i3, i4;
1052                         i1 = mol->impropers[i * 4];
1053                         i2 = mol->impropers[i * 4 + 1];
1054                         i3 = mol->impropers[i * 4 + 2];
1055                         i4 = mol->impropers[i * 4 + 3];
1056                         type1 = ATOM_AT_INDEX(mol->atoms, i1)->type;
1057                         type2 = ATOM_AT_INDEX(mol->atoms, i2)->type;
1058                         type3 = ATOM_AT_INDEX(mol->atoms, i3)->type;
1059                         type4 = ATOM_AT_INDEX(mol->atoms, i4)->type;
1060                         if (mol->par != NULL) {
1061                                 tp = ParameterLookupImproperPar(mol->par, type1, type2, type3, type4, i1, i2, i3, i4, 0);
1062                                 if (tp != NULL) {
1063                                         arena->improper_par_i[i] = (tp - mol->par->improperPars) + ATOMS_MAX_NUMBER * 2;
1064                                         continue;
1065                                 }
1066                         }
1067                         tp = ParameterLookupImproperPar(gBuiltinParameters, type1, type2, type3, type4, -1, -1, -1, -1, 0);
1068                         if (tp != NULL) {
1069                                 arena->improper_par_i[i] = (tp - gBuiltinParameters->improperPars) + ATOMS_MAX_NUMBER;
1070                                 continue;
1071                         }
1072                         /*  Record as missing  */
1073                         tp = ParameterLookupImproperPar(par, type1, type2, type3, type4, -1, -1, -1, -1, kParameterLookupMissing | kParameterLookupNoBaseAtomType);
1074                         if (tp == NULL) {
1075                                 tp = AssignArray(&par->improperPars, &par->nimproperPars, sizeof(TorsionPar), par->nimproperPars, NULL);
1076                                 tp->src = -1;
1077                                 tp->type1 = type1;
1078                                 tp->type2 = type2;
1079                                 tp->type3 = type3;
1080                                 tp->type4 = type4;
1081                         }
1082                         arena->improper_par_i[i] = (tp - par->improperPars);
1083                 }
1084                 nmissing = par->nimproperPars;
1085
1086                 /*  Copy the improper parameters  */
1087                 for (i = 0, idx = par->nimproperPars; i < mol->nimpropers; i++) {
1088                         t1 = arena->improper_par_i[i];
1089                         if (t1 < ATOMS_MAX_NUMBER)
1090                                 continue;
1091                         arena->improper_par_i[i] = idx;
1092                         for (j = i + 1; j < mol->nimpropers; j++) {
1093                                 if (arena->improper_par_i[j] == t1)
1094                                         arena->improper_par_i[j] = idx;
1095                         }
1096                         if (t1 >= ATOMS_MAX_NUMBER * 2)
1097                                 tp = mol->par->improperPars + (t1 - ATOMS_MAX_NUMBER * 2);
1098                         else tp = gBuiltinParameters->improperPars + (t1 - ATOMS_MAX_NUMBER);
1099                         AssignArray(&par->improperPars, &par->nimproperPars, sizeof(TorsionPar), idx, tp);
1100                         idx++;
1101                 }
1102         }
1103
1104         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
1105                 Int j;
1106                 char s1[8], s2[8], s3[8], s4[8];
1107                 fprintf(arena->debug_result, "\n  Improper parameters\n");
1108                 fprintf(arena->debug_result, "  No. atom1 atom2 atom3 atom4 type1 type2 type3 type4 mult  phi0      k   per\n");
1109                 for (i = 0; i < mol->nimpropers; i++) {
1110                         idx = arena->improper_par_i[i];
1111                         fprintf(arena->debug_result, "%5d %5d %5d %5d %5d ", i+1, mol->impropers[i*4]+1, mol->impropers[i*4+1]+1, mol->impropers[i*4+2]+1, mol->impropers[i*4+3]+1);
1112                         if (idx < 0) {
1113                                 fprintf(arena->debug_result, "<< missing >>\n");
1114                                 continue;
1115                         }
1116                         fprintf(arena->debug_result, "%-4s  %-4s  %-4s  %-4s  %2d ", AtomTypeDecodeToString(par->improperPars[idx].type1, s1), AtomTypeDecodeToString(par->improperPars[idx].type2, s2), AtomTypeDecodeToString(par->improperPars[idx].type3, s3), AtomTypeDecodeToString(par->improperPars[idx].type4, s4), par->improperPars[idx].mult);
1117                         for (j = 0; j < par->improperPars[idx].mult; j++) {
1118                                 fprintf(arena->debug_result, "%7.3f %7.3f %1d", par->improperPars[idx].phi0[j]*180/PI, par->improperPars[idx].k[j]*INTERNAL2KCAL, par->improperPars[idx].period[j]);
1119                         }
1120                         fprintf(arena->debug_result, "\n");
1121                 }
1122         }
1123         arena->nmissing += nmissing;
1124         return nmissing;
1125 }
1126
1127 /*  Find one fragment, starting from start_index  */
1128 static void
1129 s_find_fragment_sub(MDArena *arena, Int start_index, Int fragment_index)
1130 {
1131         Atom *atoms = arena->mol->atoms;
1132         int i, j;
1133         for (i = 0; i < atoms[start_index].nconnects; i++) {
1134                 AtomConnectEntryAtIndex(ATOM_AT_INDEX(atoms, start_index), i);
1135                 if (j >= 0 && j < arena->natoms_uniq) {
1136                         int n = arena->fragment_indices[j];
1137                         if (n < 0) {
1138                                 arena->fragment_indices[j] = fragment_index;
1139                                 s_find_fragment_sub(arena, j, fragment_index);
1140                         } else if (n != fragment_index) {
1141                                 if (arena->log_result != NULL)
1142                                         fprintf(arena->log_result, "Warning: internal inconsistency in finding fragment at atom %d\n", j + 1);
1143                         }
1144                 }
1145         }
1146 }
1147
1148 /*  Find fragments  */
1149 void
1150 md_find_fragments(MDArena *arena)
1151 {
1152         int i, idx, nuniq;
1153         if (arena->fragment_indices != NULL)
1154                 free(arena->fragment_indices);
1155         arena->fragment_indices = (Int *)malloc(sizeof(Int) * arena->natoms_uniq);
1156         if (arena->fragment_indices == NULL)
1157                 md_panic(arena, "Low memory in md_find_fragments");
1158         nuniq = arena->natoms_uniq;
1159         for (i = 0; i < nuniq; i++)
1160                 arena->fragment_indices[i] = -1;
1161         idx = 0;
1162         for (i = 0; i < nuniq; i++) {
1163                 if (arena->fragment_indices[i] >= 0)
1164                         continue;
1165                 s_find_fragment_sub(arena, i, idx);
1166                 idx++;
1167         }
1168         arena->nfragments = idx;
1169         if (arena->fragment_info != NULL)
1170                 free(arena->fragment_info);
1171         arena->fragment_info = (struct MDFragmentInfo *)calloc(sizeof(struct MDFragmentInfo), idx);
1172         if (arena->fragment_info == NULL)
1173                 md_panic(arena, "Low memory in md_find_fragments");
1174 }
1175
1176 /*  Calculate center-of-mass  */
1177 void
1178 md_center_of_mass(MDArena *arena, Vector *cp)
1179 {
1180         int i, n;
1181         Atom *ap;
1182         Double sumw;
1183         VecZero(*cp);
1184         sumw = 0.0;
1185         n = arena->mol->natoms;
1186         for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1187                 VecScaleInc(*cp, ap->r, ap->weight);
1188                 sumw += ap->weight;
1189         }
1190         VecScaleSelf(*cp, 1.0 / sumw);
1191 }
1192
1193 /*  Relocate center-of-mass  */
1194 int
1195 md_relocate_center(MDArena *arena)
1196 {
1197         int i, n;
1198         Atom *ap;
1199         Vector cm;
1200         n = arena->mol->natoms;
1201         md_center_of_mass(arena, &cm);
1202         VecDec(cm, arena->initial_center);
1203         for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1204                 VecDec(ap->r, cm);
1205         }
1206         return 0;
1207 }
1208
1209 /*  Quench translational and rotational momenta  */
1210 int
1211 md_quench_momenta(MDArena *arena)
1212 {
1213         int i, n;
1214         Atom *ap;
1215         Vector am, cm, m, v, v1, v2;
1216         Double w, sumw, kinetic;
1217         Mat33 mi;
1218
1219         if (!arena->quench_angular_momentum && !arena->quench_translational_momentum)
1220                 return 0;
1221
1222         n = arena->mol->natoms;
1223
1224         /*  Center of mass  */
1225         md_center_of_mass(arena, &cm);
1226
1227         /*  Initial kinetic energy (times 2)  */
1228         kinetic = 0.0;
1229         for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1230                 kinetic += ap->weight * VecLength2(ap->v);
1231         }
1232
1233         /*  Quench the angular momentum  */
1234         /*  Calculate the total angular momentum  */
1235
1236         if (arena->quench_angular_momentum) {
1237                 VecZero(am);
1238                 for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1239                         VecSub(v, ap->r, cm);
1240                         VecCross(v1, v, ap->v);
1241                         VecScaleInc(am, v1, ap->weight);
1242                 }
1243                 /*  Moment of inertia  */
1244                 memset(mi, 0, sizeof(mi));
1245                 for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1246                         Double x = ap->r.x - cm.x;
1247                         Double y = ap->r.y - cm.y;
1248                         Double z = ap->r.z - cm.z;
1249                         w = ap->weight;
1250                         mi[0] += w * (y * y + z * z);
1251                         mi[4] += w * (z * z + x * x);
1252                         mi[8] += w * (x * x + y * y);
1253                         mi[1] -= w * x * y;
1254                         mi[5] -= w * y * z;
1255                         mi[2] -= w * z * x;
1256                 }
1257                 mi[3] = mi[1];
1258                 mi[6] = mi[2];
1259                 mi[7] = mi[5];
1260                 /*  Calculate the angular velocity as a rigid body  */
1261                 /*  I * w = L; I, moment of inertia; w, angular velocity; L, angular momentum  */
1262                 MatrixInvert(mi, mi);
1263                 MatrixVec(&v, mi, &am);
1264                 /*  Subtract the velocity at the atom position derived from the angular velocity  */
1265                 for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1266                         VecSub(v1, ap->r, cm);
1267                         VecCross(v2, v, v1);
1268                         VecDec(ap->v, v2);
1269                 }
1270         }
1271
1272         /*  Quench the translational momentum  */
1273         if (arena->quench_translational_momentum) {
1274                 VecZero(m);
1275                 sumw = 0.0;
1276                 for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1277                         VecScaleInc(m, ap->v, ap->weight);
1278                         sumw += ap->weight;
1279                 }
1280                 VecScaleSelf(m, 1.0 / sumw);
1281                 for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++)
1282                         VecDec(ap->v, m);
1283         }
1284
1285         /*  Current kinetic energy (times 2) */
1286         w = 0.0;
1287         for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1288                 w += ap->weight * VecLength2(ap->v);
1289         }
1290         w = sqrt(kinetic / w);
1291
1292         /*  Scale the velocities to keep the kinetic energy  */
1293         for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
1294                 VecScaleSelf(ap->v, w);
1295         }
1296         return 0;
1297 }
1298
1299 /*  Initilize the runtime fields that are affected by the atom positions  */
1300 /*  (Should be called after coordinates are updated without changing structure info)  */
1301 void
1302 md_init_for_positions(MDArena *arena)
1303 {
1304         Molecule *mol;
1305         Parameter *par;
1306         int i, j, idx;
1307
1308         if (arena == NULL || (mol = arena->mol) == NULL || (par = arena->par) == NULL)
1309                 return;
1310
1311         /*  Initialize fix positions  */
1312         for (i = j = 0; i < mol->natoms; i++) {
1313                 Atom *ap = &(mol->atoms[i]);
1314                 if (ap->fix_force != 0.0) {
1315                         ap->fix_pos = ap->r;
1316                         j++;
1317                 }
1318         }
1319         if (j > 0)
1320                 md_log(arena, "Number of fixed atoms = %d\n", j);
1321 /*      if (arena->fix_atoms != NULL) {
1322                 for (i = 0; i < arena->nfix_atoms; i++) {
1323                         j = arena->fix_atoms[i].index;
1324                         if (j >= 0 && j < mol->natoms)
1325                                 arena->fix_atoms[i].pos = mol->atoms[j].r;
1326                 }
1327                 printf("Number of fixed atoms = %d\n", arena->nfix_atoms);
1328         } */
1329
1330         /*  Abnormal bonds  */
1331         if (arena->anbond_thres > 0.0) {
1332                 Vector r12;
1333                 Double r;
1334                 for (i = 0; i < mol->nbonds; i++) {
1335                         idx = arena->bond_par_i[i];
1336                         if (idx < 0)
1337                                 continue;
1338                         VecSub(r12, mol->atoms[mol->bonds[i*2]].r, mol->atoms[mol->bonds[i*2+1]].r);
1339                         r = VecLength(r12);
1340                         if (r >= (1 + arena->anbond_thres) * par->bondPars[idx].r0)
1341                                 arena->anbond_r0[i] = r - par->bondPars[idx].r0;
1342                         else
1343                                 arena->anbond_r0[i] = 0.0;
1344                 }
1345         }
1346
1347         /*  Center of mass  */
1348         md_center_of_mass(arena, &(arena->initial_center));
1349 }
1350
1351 /*  Set the alchemical flags  */
1352 /*  Independent with arena->xmol, mol  */
1353 int
1354 md_set_alchemical_flags(MDArena *arena, int nflags, const char *flags)
1355 {
1356         if (arena == NULL)
1357                 return 0;
1358
1359         if (nflags == 0 || flags == NULL) {
1360                 if (arena->alchem_flags != NULL)
1361                         free(arena->alchem_flags);
1362                 arena->alchem_flags = NULL;
1363                 arena->nalchem_flags = 0;
1364                 return 0;
1365         }
1366
1367         if (arena->alchem_flags == NULL)
1368                 arena->alchem_flags = (char *)malloc(nflags);
1369         else arena->alchem_flags = (char *)realloc(arena->alchem_flags, nflags);
1370         if (arena->alchem_flags == NULL) {
1371                 arena->nalchem_flags = 0;
1372                 return -1;
1373         }
1374         memmove(arena->alchem_flags, flags, nflags);
1375         arena->nalchem_flags = nflags;
1376         return 0;
1377 }
1378
1379 /*  Set the external forces  */
1380 /*  Independent with arena->xmol, mol  */
1381 int
1382 md_set_external_forces(MDArena *arena, int nexforces, const Vector *exforces)
1383 {
1384         if (arena == NULL)
1385                 return 0;
1386
1387         if (nexforces == 0 || exforces == NULL) {
1388                 if (arena->exforces != NULL)
1389                         free(arena->exforces);
1390                 arena->exforces = NULL;
1391                 arena->nexforces = 0;
1392                 return 0;
1393         }
1394
1395         if (arena->exforces == NULL)
1396                 arena->exforces = (Vector *)malloc(nexforces * sizeof(Vector));
1397         else arena->exforces = (Vector *)realloc(arena->exforces, nexforces * sizeof(Vector));
1398         if (arena->exforces == NULL) {
1399                 arena->nexforces = 0;
1400                 return -1;
1401         }
1402         memmove(arena->exforces, exforces, nexforces * sizeof(Vector));
1403         arena->nexforces = nexforces;
1404         return 0;
1405 }
1406
1407 const char *
1408 md_prepare(MDArena *arena, int check_only)
1409 {
1410         Int i, idx;
1411         Int t1, t2, t3;
1412         Molecule *mol;
1413
1414         if (arena->xmol->needsMDRebuild) {
1415                 /*  Set molecule again  */
1416                 md_arena_set_molecule(arena, arena->xmol);
1417                 arena->xmol->needsMDRebuild = 0;
1418         } else {
1419                 md_copy_coordinates_to_internal(arena);
1420         }
1421         mol = arena->mol;
1422
1423         arena->errmsg[0] = 0;
1424         arena->setjmp_buf = NULL;
1425
1426         /*  Table of missing parameters  */
1427 /*      Int *missing = NULL; */
1428 /*      Int nmissing = 0;  */
1429
1430         if (mol == NULL)
1431                 return "molecule is not defined";
1432         if (mol->natoms == 0 || mol->atoms == NULL)
1433                 return "molecule has no atoms";
1434
1435         /*  Open files  */
1436         if (!check_only) {
1437                 char *cwd = NULL;
1438                 const char *err = NULL;
1439
1440                 if (arena->xmol->path != NULL || mol->path != NULL) {
1441                         /*  Temporarily change to the document directory  */
1442                         char *p;
1443                         char *fname = (char *)arena->xmol->path;
1444                         if (fname == NULL)
1445                                 fname = (char *)mol->path;
1446                         fname = strdup(fname);
1447                         if ((p = strrchr(fname, '/')) != NULL
1448                                 || (p = strrchr(fname, '\\')) != NULL
1449                                 ) {
1450                                 *p = 0;
1451                                 cwd = getcwd(NULL, 0);
1452                                 chdir(fname);
1453                         }
1454                         free(fname);
1455                 }
1456                 if (arena->log_result_name != NULL && arena->log_result == NULL) {
1457                         arena->log_result = fopen(arena->log_result_name, "wb");
1458                         if (arena->log_result == NULL)
1459                                 err = "cannot create log file";
1460                 }
1461                 if (err == NULL && arena->coord_result_name != NULL && arena->coord_result == NULL) {
1462                         char molname[64];
1463                         arena->coord_result = fopen(arena->coord_result_name, "wb");
1464                         if (arena->coord_result == NULL)
1465                                 err = "cannot create coord file";
1466                         else {
1467                                 int natoms = (arena->output_expanded_atoms ? arena->mol->natoms : arena->natoms_uniq);
1468                                 MoleculeCallback_displayName(mol, molname, sizeof molname);
1469                                 fprintf(arena->coord_result, "TITLE: %s, %d atoms\n", molname, natoms);
1470                                 fflush(arena->coord_result);
1471                         }
1472                 }
1473                 if (err == NULL && arena->vel_result_name != NULL && arena->vel_result == NULL) {
1474                         arena->vel_result = fopen(arena->vel_result_name, "wb");
1475                         if (arena->vel_result == NULL)
1476                                 err = "cannot create vel file";
1477                 }
1478                 if (err == NULL && arena->force_result_name != NULL && arena->force_result == NULL) {
1479                         arena->force_result = fopen(arena->force_result_name, "wb");
1480                         if (arena->force_result == NULL)
1481                                 err = "cannot create force file";
1482                 }
1483                 if (err == NULL && arena->extend_result_name != NULL && arena->extend_result == NULL) {
1484                         arena->extend_result = fopen(arena->extend_result_name, "wb");
1485                         if (arena->extend_result == NULL)
1486                                 err = "cannot create extend file";
1487                 }
1488                 if (err == NULL && arena->debug_result_name != NULL && arena->debug_result == NULL) {
1489                         arena->debug_result = fopen(arena->debug_result_name, "wb");
1490                         if (arena->debug_result == NULL)
1491                                 err = "cannot create debug file";
1492                 }
1493                 if (cwd != NULL) {
1494                         chdir(cwd);
1495                         free(cwd);
1496                 }
1497                 if (err != NULL)
1498                         return err;
1499         }
1500
1501         /*  Count symmetry unique atoms  */
1502         /*  Atoms should be ordered so that all symmetry related atoms come after symmetry unique atoms  */
1503         for (i = t1 = 0, t2 = -1; i < mol->natoms; i++) {
1504                 Atom *ap = mol->atoms + i;
1505                 Symop symop = ap->symop;
1506                 if (SYMOP_ALIVE(symop)) {
1507                         if (t2 == -1)
1508                                 t2 = i;  /*  The index of the first non-unique atom  */
1509                 } else {
1510                         t1++;        /*  The number of unique atoms  */
1511                 }
1512         }
1513         if (t2 >= 0 && t1 != t2)
1514                 return "all symmetry related atoms should be after symmetry unique atoms";
1515 /*      if (t1 > mol->natoms && mol->box == NULL)
1516                 return "symmetry operation is used but no periodic box is defined"; */
1517
1518         md_set_cell(arena);
1519
1520         arena->natoms_uniq = t1;
1521
1522         if (arena->debug_result != NULL) {
1523                 time_t loc_time;
1524                 time(&loc_time);
1525                 fprintf(arena->debug_result, "---- LWMD Started at %s", ctime(&loc_time));
1526         }
1527
1528         /*  Recalc angle/dihedral/improper tables from the bond table and parameters  */
1529 /*      s_rebuild_angles(arena);
1530         s_rebuild_dihedrals(arena);
1531         s_rebuild_impropers(arena); */
1532
1533         /*  Statistics of the molecule  */
1534         md_log(arena,
1535                    "Number of bonds = %d\n"
1536                    "Number of angles = %d\n"
1537                    "Number of dihedrals = %d\n"
1538                    "Number of impropers = %d\n", mol->nbonds, mol->nangles, mol->ndihedrals, mol->nimpropers);
1539
1540         t1 = t2 = t3 = 0;
1541         for (i = 0; i < mol->natoms; i++) {
1542                 if (mol->atoms[i].fix_force > 0)
1543                         t1++;
1544                 if (mol->atoms[i].fix_force < 0)
1545                         t2++;
1546                 if (mol->atoms[i].mm_exclude)
1547                         t3++;
1548         }
1549         if (t1 > 0)
1550                 md_log(arena, "Number of constrained atoms = %d\n", t1);
1551         if (t2 > 0)
1552                 md_log(arena, "Number of fixed atoms = %d\n", t2);
1553         if (t3 > 0)
1554                 md_log(arena, "Number of excluded atoms = %d\n", t3);
1555
1556         if (arena->natoms_uniq < mol->natoms) {
1557                 md_log(arena, "Number of symmetry-unique atoms = %d\n", arena->natoms_uniq);
1558                 t2 = 0;
1559                 for (i = 0; i < arena->natoms_uniq; i++) {
1560                         if (mol->atoms[i].fix_force < 0)
1561                                 t2++;
1562                 }
1563         }
1564
1565         arena->degree_of_freedom = 3 * (arena->natoms_uniq - t2 - t3);
1566         md_log(arena, "Degree of freedom = %d\n", arena->degree_of_freedom);
1567
1568         /*  Build local cache of the used parameters  */
1569         if (arena->par != NULL)
1570                 ParameterRelease(arena->par);
1571         arena->par = ParameterNew();
1572         arena->nmissing = arena->nsuspicious = 0;
1573         s_find_vdw_parameters(arena);
1574         s_find_bond_parameters(arena);
1575         s_find_angle_parameters(arena);
1576         s_find_dihedral_parameters(arena);
1577         s_find_improper_parameters(arena);
1578
1579         if (arena->nmissing > 0) {
1580         /*      for (i = 0; i < nmissing; i++) {
1581                         char s1[6], s2[6], s3[6], s4[6];
1582                         type1 = missing[i * 5];
1583                         AtomTypeDecodeToString(missing[i * 5 + 1], s1);
1584                         AtomTypeDecodeToString(missing[i * 5 + 2], s2);
1585                         AtomTypeDecodeToString(missing[i * 5 + 3], s3);
1586                         AtomTypeDecodeToString(missing[i * 5 + 4], s4);
1587                         if (type1 == 0)
1588                                 md_warning(arena, "Missing vdw parameter for %s\n", s1, s2);
1589                         else if (type1 == 1)
1590                                 md_warning(arena, "Missing bond parameter for %s-%s\n", s1, s2);
1591                         else if (type1 == 2)
1592                                 md_warning(arena, "Missing angle parameter for %s-%s-%s\n", s1, s2, s3);
1593                         else if (type1 == 3)
1594                                 md_warning(arena, "Missing dihedral parameter for %s-%s-%s-%s\n", s1, s2, s3, s4);
1595                         else if (type1 == 4)
1596                                 md_warning(arena, "Missing improper parameter for %s-%s-%s-%s\n", s1, s2, s3, s4);
1597                 }
1598                 free(missing); */
1599                 return "some parameters are missing";
1600         }
1601
1602         /*  Build the exclusion table  */
1603         s_make_exclusion_list(arena);
1604         if (arena->debug_result != NULL && arena->debug_output_level > 0) {
1605                 fprintf(arena->debug_result, "\n  MDExclusion table for each atom\n");
1606                 fprintf(arena->debug_result, "  No.  {self and special} {1-2} {1-3} {1-4}\n");
1607                 for (i = 0; i < mol->natoms; i++) {
1608                         fprintf(arena->debug_result, "%5d ", i+1);
1609                         for (idx = arena->exinfo[i].index0; idx <= arena->exinfo[i+1].index0; idx++) {
1610                                 int n = 0;
1611                                 if (idx == arena->exinfo[i].index0)
1612                                         n += fprintf(arena->debug_result, "{");
1613                                 if (idx == arena->exinfo[i].index1)
1614                                         n += fprintf(arena->debug_result, "} {");
1615                                 if (idx == arena->exinfo[i].index2)
1616                                         n += fprintf(arena->debug_result, "} {");
1617                                 if (idx == arena->exinfo[i].index3)
1618                                         n += fprintf(arena->debug_result, "} {");
1619                                 if (idx < arena->exinfo[i+1].index0) {
1620                                         if (n == 0)
1621                                                 fprintf(arena->debug_result, " ");
1622                                         fprintf(arena->debug_result, "%d", arena->exlist[idx]+1);
1623                                 }
1624                         }
1625                         fprintf(arena->debug_result, "}\n");
1626                 }
1627         }
1628
1629         /*  Parameter checking only  */
1630         if (check_only) {
1631                 arena->is_initialized = 1;  /*  Only static fields are ready  */
1632                 arena->mol->needsMDRebuild = 0;
1633                 return NULL;
1634         }
1635
1636         /*  Allocate storage for Verlet list  */
1637         arena->max_nverlets = mol->natoms;
1638         if (arena->verlets != NULL)
1639                 free(arena->verlets);
1640         arena->verlets = (MDVerlet *)calloc(sizeof(MDVerlet), arena->max_nverlets);
1641         arena->nverlets = 0;
1642         if (arena->verlets_dr != NULL)
1643                 free(arena->verlets_dr);
1644         arena->verlets_dr = (Vector *)calloc(sizeof(Vector), mol->natoms);
1645         if (arena->verlets == NULL || arena->verlets_dr == NULL)
1646                 md_panic(arena, ERROR_out_of_memory);
1647         arena->last_verlet_step = -1;
1648         if (arena->verlet_i != NULL)
1649                 free(arena->verlet_i);
1650         arena->verlet_i = (Int *)calloc(sizeof(Int), mol->natoms + 1);
1651
1652         /*  Allocate storage for partial energy/force  */
1653         if (arena->energies != NULL)
1654                 free(arena->energies);
1655         arena->energies = (Double *)calloc(sizeof(Double), kEndIndex);
1656         if (arena->forces != NULL)
1657                 free(arena->forces);
1658         arena->forces = (Vector *)calloc(sizeof(Vector), kKineticIndex * mol->natoms);
1659         if (arena->energies == NULL || arena->forces == NULL)
1660                 md_panic(arena, ERROR_out_of_memory);
1661
1662         /*  Initialize storage for pair interaction calculations */
1663         /*  (The storage will be allocated when necessary)  */
1664 /*      if (arena->group_flags != NULL) {
1665                 free(arena->group_flags);
1666                 arena->group_flags = NULL;
1667         } */
1668         if (arena->pair_forces != NULL) {
1669                 free(arena->pair_forces);
1670                 arena->pair_forces = NULL;
1671         }
1672
1673         /*  Allocate ring buffer   */
1674         if (arena->ring != NULL) {
1675                 free(arena->ring->buf);
1676                 free(arena->ring);
1677         }
1678         arena->ring = (MDRing *)calloc(sizeof(MDRing), 1);
1679         if (arena->ring == NULL)
1680                 md_panic(arena, ERROR_out_of_memory);
1681         arena->ring->size = mol->natoms;
1682         if (arena->pressure != NULL && arena->pressure->disabled == 0)
1683                 arena->ring->size = mol->natoms + 4;
1684 #if MINIMIZE_CELL
1685         if (arena->minimize_cell && arena->mol->cell != NULL)
1686                 arena->ring->size = mol->natoms + 4;
1687 #endif
1688         arena->ring->nframes = 2000 / arena->ring->size;
1689         if (arena->ring->nframes < 2)
1690                 arena->ring->nframes = 2;
1691         arena->ring->buf = (Vector *)calloc(sizeof(Vector), arena->ring->size * arena->ring->nframes);
1692         if (arena->ring->buf == NULL)
1693                 md_panic(arena, ERROR_out_of_memory);
1694         arena->ring->next = 0;
1695         arena->ring->count = 0;
1696
1697         /*  Initialize temperature statistics  */
1698         arena->sum_temperature = 0.0;
1699         arena->nsum_temperature = 0;
1700
1701         /*  Initialize unique bond/angle/dihedral/improper table  */
1702 /*      s_find_symmetry_unique(arena); */
1703
1704         /*  Initialize the position-dependent fields  */
1705         md_init_for_positions(arena);
1706
1707         /*  Clear the snapshot storage  */
1708 /*      if (arena->snapshots != NULL) {
1709                 for (i = 0; i < arena->nsnapshots; i++) {
1710                         if (arena->snapshots[i] != NULL)
1711                                 free(arena->snapshots[i]);
1712                 }
1713                 free(arena->snapshots);
1714                 arena->snapshots = NULL;
1715         }
1716         arena->nsnapshots = 0;
1717 */
1718
1719         /*  Random number seed  */
1720         {
1721                 unsigned int seed = arena->random_seed;
1722                 if (seed == 0)
1723                         seed = (unsigned int)time(NULL);
1724                 md_srand(seed);
1725                 md_log(arena, "Random number seed = %u\n", seed);
1726         }
1727
1728         /*  Clear the surface area record (will be reallocated within calc_surface_force() if necessary) */
1729         if (arena->sp_arena != NULL) {
1730                 clear_sp_arena(arena);
1731         }
1732 /*      if (arena->sp2_arena != NULL) {
1733                 clear_sp2_arena(arena);
1734         }
1735  */
1736
1737         /*  Graphite potential  */
1738         if (arena->use_graphite) {
1739                 if (arena->graphite == NULL)
1740                         arena->graphite = graphite_new();
1741                 graphite_set_needs_update(arena->graphite, 1);
1742         }
1743
1744         if (arena->velocities_read == 0)
1745                 md_init_velocities(arena);
1746
1747         /*  Fragment analysis  */
1748         md_find_fragments(arena);
1749         md_log(arena, "%d fragments found\n", arena->nfragments);
1750
1751         /*  Pressure control statistics  */
1752         if (arena->pressure != NULL) {
1753                 if (mol->cell == NULL)
1754                         return "pressure control is requested but no unit cell is defined";
1755                 pressure_prepare(arena);
1756         }
1757
1758         arena->is_initialized = 2;   /*  Runtime fields are ready  */
1759         arena->mol->needsMDRebuild = 0;
1760         arena->request_abort = 0;
1761
1762         return NULL;
1763 }
1764
1765 #if __WXMSW__
1766 #define random rand
1767 #define srandom srand
1768 #define kRandMax ((double)RAND_MAX)
1769 #else
1770 #define kRandMax 2147483648.0
1771 #endif
1772
1773 /*  A uniform random number in [0,1]  */
1774 Double
1775 md_rand(void)
1776 {
1777         return (double)random() / kRandMax;
1778 }
1779
1780 /*  A random number with gaussian distribution  */
1781 Double
1782 md_gaussian_rand(void)
1783 {
1784         static int waiting = 0;
1785         static Double next_value = 0;
1786         Double f, r, v1, v2;
1787         if (waiting) {
1788                 waiting = 0;
1789                 return next_value;
1790         }
1791         r = 2.0;
1792         while (r >= 1.0 || r < 1.523e-8) {
1793                 v1 = 2.0 * md_rand() - 1.0;
1794                 v2 = 2.0 * md_rand() - 1.0;
1795                 r = v1 * v1 + v2 * v2;
1796         }
1797         f = sqrt(-2.0 * log(r) / r);
1798         waiting = 1;
1799         next_value = v1 * f;
1800         return v2 * f;
1801 }
1802
1803 /*  Seed the random number  */
1804 void
1805 md_srand(unsigned int seed)
1806 {
1807         srandom(seed);
1808 }
1809
1810 /*  Scale velocities to match the given temperature  */
1811 void
1812 md_scale_velocities(MDArena *arena)
1813 {
1814         int i;
1815         Atom *ap;
1816         Double ttemp, scale;
1817         if (arena->nsum_temperature == 0) {
1818                 Double kinetic = 0.0;
1819                 for (i = 0, ap = arena->mol->atoms; i < arena->natoms_uniq; i++, ap++) {
1820                         if (ap->fix_force < 0 || ap->mm_exclude)
1821                                 continue;
1822                         kinetic += ap->weight * VecLength2(ap->v);
1823                 }
1824                 kinetic *= 0.5;
1825                 ttemp = 2.0 * kinetic / (arena->degree_of_freedom * BOLTZMANN);
1826         } else {
1827                 ttemp = arena->sum_temperature / arena->nsum_temperature;
1828         }
1829         scale = sqrt(arena->temperature / ttemp);
1830         for (i = 0, ap = arena->mol->atoms; i < arena->mol->natoms; i++, ap++) {
1831                 if (ap->fix_force < 0)
1832                         continue;
1833                 VecScaleSelf(ap->v, scale);
1834         }
1835         arena->sum_temperature = 0;
1836         arena->nsum_temperature = 0;
1837 }
1838
1839 /*  Give random velocities that matches the given temperature  */
1840 void
1841 md_init_velocities(MDArena *arena)
1842 {
1843         int i, n;
1844         Double w;
1845 /*      Double temp = arena->temperature; */
1846         Atom *ap = arena->mol->atoms;
1847         n = arena->mol->natoms;
1848         for (i = 0; i < n; i++, ap++) {
1849                 if (ap->fix_force < 0 || fabs(ap->weight) < 1e-6 || ap->mm_exclude) {
1850                         ap->v.x = ap->v.y = ap->v.z = 0;
1851                 } else {
1852                         w = sqrt(arena->temperature * BOLTZMANN / ap->weight);
1853                         ap->v.x = w * md_gaussian_rand();
1854                         ap->v.y = w * md_gaussian_rand();
1855                         ap->v.z = w * md_gaussian_rand();
1856                 /*      ap->v.x = md_rand() - 0.5;
1857                         ap->v.y = md_rand() - 0.5;
1858                         ap->v.z = md_rand() - 0.5; */
1859                 }
1860         }
1861         arena->sum_temperature = 0;
1862         arena->nsum_temperature = 0;
1863
1864         /*  Quench the total momentum  */
1865         md_quench_momenta(arena);
1866
1867         /*  Adjust the temperature  */
1868         md_scale_velocities(arena);
1869
1870         arena->velocities_read = 0;
1871 }
1872
1873 void
1874 md_bootstrap(MDArena *arena)
1875 {
1876         /*  Set the initial velocities and calculate the current force  */
1877         md_init_velocities(arena);
1878         calc_force(arena);
1879 }
1880
1881 static int
1882 s_md_output_extend_info(MDArena *arena, FILE *fp)
1883 {
1884         Vector *ap, *bp, *cp, *op;
1885         static Vector zero = {0.0, 0.0, 0.0};
1886         if (arena->coord_result_frame % 10 == 0 || ftello(fp) == 0) {
1887                 fprintf(fp, "EXLABEL: %7s %7s %11s ", "FRAME", "STEP", "POT_ENERGY");
1888                 fprintf(fp, "%7s %7s %7s ", "CELL_AX", "CELL_AY", "CELL_AZ");
1889                 fprintf(fp, "%7s %7s %7s ", "CELL_BX", "CELL_BY", "CELL_BZ");
1890                 fprintf(fp, "%7s %7s %7s ", "CELL_CX", "CELL_CY", "CELL_CZ");
1891                 fprintf(fp, "%7s %7s %7s\n", "CELL_OX", "CELL_OY", "CELL_OZ");
1892         }
1893
1894         fprintf(fp, "EXINFO:  %7d %7d %11.5f ", arena->coord_result_frame, arena->step, arena->total_energy * INTERNAL2KCAL);
1895         if (arena->mol->cell != NULL) {
1896                 ap = &(arena->mol->cell->axes[0]);
1897                 bp = &(arena->mol->cell->axes[1]);
1898                 cp = &(arena->mol->cell->axes[2]);
1899                 op = &(arena->mol->cell->origin);
1900         } else {
1901                 ap = bp = cp = op = &zero;
1902         }
1903         fprintf(fp, "%7.3f %7.3f %7.3f ", ap->x, ap->y, ap->z);
1904         fprintf(fp, "%7.3f %7.3f %7.3f ", bp->x, bp->y, bp->z);
1905         fprintf(fp, "%7.3f %7.3f %7.3f ", cp->x, cp->y, cp->z);
1906         fprintf(fp, "%7.3f %7.3f %7.3f\n", op->x, op->y, op->z);
1907
1908         fflush(fp);
1909
1910         return 0;
1911 }
1912
1913 int
1914 md_output_results(MDArena *arena)
1915 {
1916         int i, j, natoms;
1917         Atom *ap;
1918         Vector *av, *bv, *cv;
1919
1920         natoms = (arena->output_expanded_atoms ? arena->mol->natoms : arena->natoms_uniq);
1921
1922 #if 0
1923         if (arena->coord_result == NULL && arena->coord_result_name != NULL) {
1924                 arena->coord_result = fopen(arena->coord_result_name, "wb");
1925                 if (arena->coord_result == NULL) {
1926                         snprintf(sErrBuf, sizeof sErrBuf, "Cannot open coordinate result file %s", arena->coord_result_name);
1927                         arena->errmsg = sErrBuf;
1928                         return 1;
1929                 }
1930                 fprintf(arena->coord_result, "TITLE: coordinates for %d atoms\n", natoms);
1931                 arena->coord_result_frame = 0;
1932         }
1933         if (arena->vel_result == NULL && arena->vel_result_name != NULL) {
1934                 arena->vel_result = fopen(arena->vel_result_name, "wb");
1935                 if (arena->vel_result == NULL) {
1936                         snprintf(sErrBuf, sizeof sErrBuf, "Cannot open velocity result file %s", arena->vel_result_name);
1937                         arena->errmsg = sErrBuf;
1938                         return 1;
1939                 }
1940                 fprintf(arena->vel_result, "TITLE: velocities for %d atoms\n", natoms);
1941         }
1942         if (arena->force_result == NULL && arena->force_result_name != NULL) {
1943                 arena->force_result = fopen(arena->force_result_name, "wb");
1944                 if (arena->force_result == NULL) {
1945                         snprintf(sErrBuf, sizeof sErrBuf, "Cannot open force result file %s", arena->force_result_name);
1946                         arena->errmsg = sErrBuf;
1947                         return 1;
1948                 }
1949                 fprintf(arena->force_result, "TITLE: force for %d atoms\n", natoms);
1950         }
1951         if (arena->extend_result == NULL && arena->extend_result_name != NULL) {
1952                 arena->extend_result = fopen(arena->extend_result_name, "wb");
1953                 if (arena->extend_result == NULL) {
1954                         snprintf(sErrBuf, sizeof sErrBuf, "Cannot open extend info file %s", arena->extend_result_name);
1955                         arena->errmsg = sErrBuf;
1956                         return 1;
1957                 }
1958                 fprintf(arena->extend_result, "# Frame step energy ax ay az bx by bz cx cy cz ox oy oz\n");
1959         }
1960 #endif
1961
1962         ap = arena->mol->atoms;
1963         if (arena->wrap_coordinates && arena->mol->cell != NULL) {
1964                 av = &(arena->mol->cell->axes[0]);
1965                 bv = &(arena->mol->cell->axes[1]);
1966                 cv = &(arena->mol->cell->axes[2]);
1967         } else {
1968                 av = bv = cv = NULL;
1969         }
1970         if (arena->coord_result != NULL) {
1971                 if (arena->wrap_coordinates)
1972                         md_wrap_coordinates(arena);
1973                 for (i = j = 0; i < natoms; i++, j = (j + 3) % 10) {
1974                         Vector r = ap[i].r;
1975                         if (av != NULL) {
1976                                 VecScaleInc(r, *av, ap[i].wrap_dx);
1977                                 VecScaleInc(r, *bv, ap[i].wrap_dy);
1978                                 VecScaleInc(r, *cv, ap[i].wrap_dz);
1979                         }
1980                         fprintf(arena->coord_result, "%8.3f%s%8.3f%s%8.3f%s",
1981                                 r.x, (j == 9 ? "\n" : ""),
1982                                 r.y, (j == 8 ? "\n" : ""),
1983                                 r.z, (j == 7 ? "\n" : ""));
1984                 }
1985                 if (j != 0)
1986                         fprintf(arena->coord_result, "\n");
1987                 if (arena->periodic_a || arena->periodic_b || arena->periodic_c) {
1988                         fprintf(arena->coord_result, "%8.3f%8.3f%8.3f\n", arena->mol->cell->cell[0], arena->mol->cell->cell[1], arena->mol->cell->cell[2]);
1989                 }
1990                 fflush(arena->coord_result);
1991         }
1992         if (arena->vel_result != NULL) {
1993                 for (i = j = 0; i < natoms; i++, j = (j + 3) % 10) {
1994                         fprintf(arena->vel_result, " %7.3f%s %7.3f%s %7.3f%s",
1995                                 ap[i].v.x*1000, (j == 9 ? "\n" : ""),
1996                                 ap[i].v.y*1000, (j == 8 ? "\n" : ""),
1997                                 ap[i].v.z*1000, (j == 7 ? "\n" : ""));
1998                 }
1999                 if (j != 0)
2000                         fprintf(arena->vel_result, "\n");
2001                 fflush(arena->vel_result);
2002         }
2003         if (arena->force_result != NULL) {
2004                 for (i = j = 0; i < natoms; i++, j = (j + 3) % 10) {
2005                         fprintf(arena->force_result, " %7.3f%s %7.3f%s %7.3f%s",
2006                                 ap[i].f.x, (j == 9 ? "\n" : ""),
2007                                 ap[i].f.y, (j == 8 ? "\n" : ""),
2008                                 ap[i].f.z, (j == 7 ? "\n" : ""));
2009                 }
2010                 if (j != 0)
2011                         fprintf(arena->force_result, "\n");
2012                 fflush(arena->force_result);
2013         }
2014         if (arena->extend_result != NULL) {
2015                 s_md_output_extend_info(arena, arena->extend_result);
2016         } else if (arena->log_result != NULL) {
2017                 s_md_output_extend_info(arena, arena->log_result);
2018         }
2019         arena->coord_result_frame++;
2020         return 0;
2021 }
2022
2023 void
2024 md_output_energies(MDArena *arena)
2025 {
2026         int i;
2027         int periodic = (arena->mol->cell != NULL) && (arena->periodic_a || arena->periodic_b || arena->periodic_c);
2028 /*      if (arena->log_result == NULL)
2029                 return; */
2030         if ((arena->step / arena->energy_output_freq) % 10 == 0) {
2031                 md_log(arena, "ELABEL:  %11s %11s %11s %11s %11s %11s", "STEP", "TOTAL_POT", "BOND", "ANGLE", "DIHEDRAL", "IMPROPER");
2032                 md_log(arena, " %11s %11s %11s %11s %11s %11s %11s %11s", "VDW", "ELECT", "AUX", "SURFACE", "KINETIC", "NET", "TEMP", "TEMP_AVG");
2033                 if (periodic)
2034                         md_log(arena, " %11s", "VOLUME");
2035                 if (arena->nalchem_flags > 0)
2036                         md_log(arena, " %11s %11s", "LAMBDA", "DEDL");
2037                 md_log(arena, "\n");
2038         }
2039         md_log(arena, "ENERGY:  %11d %11.5f", arena->step, arena->total_energy * INTERNAL2KCAL);
2040         for (i = 0; i < kEndIndex; i++)
2041                 md_log(arena, " %11.5f", arena->energies[i] * INTERNAL2KCAL);
2042         md_log(arena, " %11.5f", (arena->energies[kKineticIndex] + arena->total_energy) * INTERNAL2KCAL);
2043         md_log(arena, " %11.5f", arena->transient_temperature);
2044         md_log(arena, " %11.5f", (arena->nsum_temperature > 0 ? arena->sum_temperature / arena->nsum_temperature : 0.0));
2045         if (periodic) {
2046                 Vector v, *av, *bv, *cv;
2047                 av = &(arena->mol->cell->axes[0]);
2048                 bv = &(arena->mol->cell->axes[1]);
2049                 cv = &(arena->mol->cell->axes[2]);
2050                 VecCross(v, *av, *bv);
2051                 md_log(arena, " %11.5f", VecDot(v, *cv));
2052         }
2053         if (arena->nalchem_flags > 0) {
2054                 md_log(arena, " %11.5f %11.5f", arena->alchem_lambda, arena->alchem_energy * INTERNAL2KCAL / arena->alchem_dlambda);
2055         }
2056         md_log(arena, "\n");
2057         md_log(arena, NULL);
2058 }
2059
2060 void
2061 md_update_velocities(MDArena *arena)
2062 {
2063         int i, natoms;
2064         Byte use_sym;
2065         Atom *ap;
2066         Double w, wt;
2067 /*      Double wsum; */
2068 /*      Vector m1, m2; */
2069         Double halftimestep = arena->timestep * 0.5;
2070         ap = arena->mol->atoms;
2071         natoms = arena->mol->natoms;
2072         use_sym = (arena->mol->nsyms > 0);
2073         Double limit = arena->velocity_limit;
2074         limit *= limit;
2075 /*      VecZero(m1);
2076         VecZero(m2);
2077         wsum = 0.0; */
2078         for (i = 0; i < natoms; i++, ap++) {
2079                 if (use_sym && ap->symop.alive)
2080                         continue;
2081                 if (ap->fix_force < 0)
2082                         continue;
2083                 if (ap->mm_exclude)
2084                         continue;
2085                 wt = ap->weight;
2086                 if (fabs(wt) < 1e-6)
2087                         continue;
2088                 w = halftimestep / wt;
2089         /*      VecScaleInc(m1, ap->v, wt); */
2090                 VecScaleInc(ap->v, ap->f, w);
2091                 w = VecLength2(ap->v);
2092                 if (w > limit || !isfinite(w))  /*  Is isfinite() available in other platforms?? */
2093                         md_panic(arena, "the velocity for atom %d exceeded limit at step %d", i+1, arena->step);
2094         /*      VecScaleInc(m2, ap->v, wt); */
2095         /*      wsum += wt; */
2096         }
2097 /*
2098         VecDec(m2, m1);
2099         w = 1.0 / wsum;
2100         VecScaleSelf(m2, w);
2101         for (i = 0, ap = arena->mol->atoms; i < natoms; i++, ap++) {
2102                 if (use_sym && ap->symop.alive)
2103                         continue;
2104                 if (ap->fix_force < 0)
2105                         continue;
2106                 VecDec(ap->v, m2);
2107         }
2108 */
2109         arena->velocities_read = 0;
2110         /*  For debug  */
2111 /*      VecZero(m2);
2112         for (i = 0, ap = arena->mol->atoms; i < natoms; i++, ap++) {
2113                 if (use_sym && ap->sym_op > 0)
2114                         continue;
2115                 wt = ap->weight;
2116                 VecScaleInc(m2, ap->v, wt);
2117         } */
2118 }
2119
2120 void
2121 md_update_positions(MDArena *arena)
2122 {
2123         int i, natoms;
2124         Atom *ap;
2125         Double timestep = arena->timestep;
2126         Vector *vdr = arena->verlets_dr;
2127         Vector dr;
2128 /*      Double w, limit;
2129         Double kinetic, kinetic_uniq; */
2130         Byte use_sym;
2131         ap = arena->mol->atoms;
2132         natoms = arena->mol->natoms;
2133 /*      limit = arena->velocity_limit;
2134         limit *= limit; */
2135         use_sym = (arena->mol->nsyms > 0);
2136 /*      kinetic = kinetic_uniq = 0.0; */
2137         for (i = 0; i < natoms; i++, ap++, vdr++) {
2138                 if (use_sym && ap->symop.alive)
2139                         continue;
2140                 if (ap->fix_force < 0)
2141                         continue;
2142                 if (ap->mm_exclude)
2143                         continue;
2144                 VecScale(dr, ap->v, timestep);
2145                 VecInc(ap->r, dr);
2146                 VecInc(*vdr, dr);
2147         }
2148
2149         /*  Update the abnormal bond parameters  */
2150         if (arena->anbond_thres > 0.0) {
2151                 Double *fp = arena->anbond_r0;
2152                 for (i = 0; i < arena->mol->nbonds; i++) {
2153                         if (fp[i] <= 0.0)
2154                                 continue;
2155                         fp[i] -= arena->anbond_anneal_rate;
2156                         if (fp[i] < 0.0)
2157                                 fp[i] = 0.0;
2158                 }
2159         }
2160
2161         /*  Relocate the center of mass (if necessary)  */
2162         if (arena->relocate_center) {
2163                 md_relocate_center(arena);
2164         }
2165 }
2166
2167 void
2168 md_calc_kinetic_energy(MDArena *arena)
2169 {
2170         int i, n, nuniq;
2171         Double w, kinetic, kinetic_uniq;
2172         Atom *ap;
2173         n = arena->mol->natoms;
2174         nuniq = arena->natoms_uniq;
2175         kinetic = kinetic_uniq = 0;
2176         for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
2177                 if (ap->fix_force < 0)
2178                         continue;
2179                 w = VecLength2(ap->v) * ap->weight;
2180                 kinetic += w;
2181                 if (i < nuniq)
2182                         kinetic_uniq += w;
2183         }
2184         arena->energies[kKineticIndex] = kinetic * 0.5;
2185         arena->transient_temperature = kinetic_uniq / (arena->degree_of_freedom * BOLTZMANN);
2186         arena->sum_temperature += arena->transient_temperature;
2187         arena->nsum_temperature++;
2188         arena->average_temperature = arena->sum_temperature / arena->nsum_temperature;
2189 }
2190
2191 /*  Andersen thermostat  */
2192 void
2193 md_andersen_thermostat(MDArena *arena)
2194 {
2195         int n = arena->natoms_uniq;
2196         int i;
2197         Double w, de;
2198         Atom *ap;
2199         Double q = arena->andersen_thermo_coupling;
2200         de = 0.0;
2201         for (i = 0, ap = arena->mol->atoms; i < n; i++, ap++) {
2202                 if (ap->fix_force < 0 || fabs(ap->weight) < 1e-6) {
2203                         ap->v.x = ap->v.y = ap->v.z = 0;
2204                 } else if (md_rand() < q) {
2205                         de -= VecLength2(ap->v) * ap->weight;
2206                         w = sqrt(arena->temperature * BOLTZMANN / ap->weight);
2207                         ap->v.x = w * md_gaussian_rand();
2208                         ap->v.y = w * md_gaussian_rand();
2209                         ap->v.z = w * md_gaussian_rand();
2210                         de += VecLength2(ap->v) * ap->weight;
2211                 }
2212         }
2213         de *= 0.5;
2214         arena->energies[kKineticIndex] += de;
2215 }
2216
2217 void
2218 md_transform_vec_by_symmetry(MDArena *arena, Vector *dst, const Vector *src, Symop rec, int no_translation)
2219 {
2220         Transform temp;
2221         if (!rec.alive)
2222                 return;
2223         if (rec.sym > 0 && rec.sym <= arena->mol->nsyms) {
2224                 memmove(temp, arena->cellsyms[rec.sym], sizeof(temp));
2225                 if (no_translation)
2226                         temp[9] = temp[10] = temp[11] = 0.0;
2227                 TransformVec(dst, temp, src);
2228         } else *dst = *src;
2229         if (!no_translation) {
2230                 Vector *vp = arena->mol->cell->axes;
2231                 VecScaleInc(*dst, vp[0], rec.dx);
2232                 VecScaleInc(*dst, vp[1], rec.dy);
2233                 VecScaleInc(*dst, vp[2], rec.dz);
2234         }
2235 }
2236
2237 void
2238 md_wrap_coordinates(MDArena *arena)
2239 {
2240         /*  Calculate the offset for each fragment to wrap into the unit cell (0,0,0)-(1,1,1) */
2241         int i, n;
2242         int last_n;
2243         Symop last_symop;
2244         Molecule *mol = arena->mol;
2245         Atom *ap;
2246
2247         if (mol == NULL || mol->natoms == 0)
2248                 return;
2249
2250         if (arena->nfragments == 0)
2251                 md_find_fragments(arena);
2252
2253         /*  Calculate the center of mass for each fragment  */
2254         for (n = 0; n < arena->nfragments; n++) {
2255                 VecZero(arena->fragment_info[n].pos);
2256                 arena->fragment_info[n].mass = 0.0;
2257         }
2258         for (i = 0, ap = mol->atoms; i < arena->natoms_uniq; i++, ap++) {
2259                 n = arena->fragment_indices[i];
2260                 VecScaleInc(arena->fragment_info[n].pos, ap->r, ap->weight);
2261                 arena->fragment_info[n].mass += ap->weight;
2262         }
2263         for (n = 0; n < arena->nfragments; n++) {
2264                 VecScaleSelf(arena->fragment_info[n].pos, 1.0/arena->fragment_info[n].mass);
2265         }
2266
2267         /*  Calculate the offset  */
2268         /*  A trick: atoms belonging to one fragment are almost always located in a consequent
2269                 positions. So we cache the 'last' information to avoid duplicate calculation
2270                 efficiently.  */
2271         last_n = -1;
2272         last_symop.alive = 0;
2273         for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap++) {
2274                 if (ap->symop.alive)
2275                         n = arena->fragment_indices[ap->symbase];
2276                 else n = arena->fragment_indices[i];
2277                 if (n == last_n && ap->symop.alive == last_symop.alive
2278                 && (!last_symop.alive || (ap->symop.dx == last_symop.dx && ap->symop.dy == last_symop.dy && ap->symop.dz == last_symop.dz && ap->symop.sym == last_symop.sym))) {
2279                         /*  Belongs to the 'last' fragment  */
2280                         ap->wrap_dx = ap[-1].wrap_dx;
2281                         ap->wrap_dy = ap[-1].wrap_dy;
2282                         ap->wrap_dz = ap[-1].wrap_dz;
2283                 } else {
2284                         /*  Calculate the offset from the position of the center of mass */
2285                         Vector r = arena->fragment_info[n].pos;
2286                         TransformVec(&r, arena->mol->cell->rtr, &r);
2287                         if (ap->symop.alive && ap->symop.sym > 0 && ap->symop.sym <= arena->mol->nsyms) {
2288                                 /*  The translational components of symop are not included  */
2289                                 TransformVec(&r, arena->mol->syms[ap->symop.sym - 1], &r);
2290                         }
2291                         ap->wrap_dx = (arena->periodic_a ? -floor(r.x) : 0);
2292                         ap->wrap_dy = (arena->periodic_b ? -floor(r.y) : 0);
2293                         ap->wrap_dz = (arena->periodic_c ? -floor(r.z) : 0);
2294                 }
2295                 last_n = n;
2296                 last_symop = ap->symop;
2297         }
2298 }
2299
2300 void
2301 md_amend_by_symmetry(MDArena *arena)
2302 {
2303         int i, natoms;
2304         Atom *ap, *ap0;
2305         if (arena->mol->nsyms == 0)
2306                 return;
2307         ap = ap0 = arena->mol->atoms;
2308         natoms = arena->mol->natoms;
2309         for (i = 0; i < natoms; i++, ap++) {
2310                 Symop symop;
2311                 if (!ap->symop.alive)
2312                         continue;
2313                 symop = ap->symop;
2314                 ap0 = arena->mol->atoms + ap->symbase;
2315                 md_transform_vec_by_symmetry(arena, &(ap->r), &(ap0->r), symop, 0);
2316                 md_transform_vec_by_symmetry(arena, &(ap->f), &(ap0->f), symop, 1);
2317                 md_transform_vec_by_symmetry(arena, &(ap->v), &(ap0->v), symop, 1);
2318         }
2319
2320 }
2321
2322 void
2323 md_snapshot(MDArena *arena, int idx)
2324 {
2325         Atom *ap;
2326         Vector *vp;
2327         MDSnapshot **spp;
2328         size_t size;
2329         int i, natoms;
2330         if (idx < 0)
2331                 return;
2332         spp = (MDSnapshot **)AssignArray(&arena->snapshots, &arena->nsnapshots, sizeof(*spp), idx, NULL);
2333         if (spp == NULL)
2334                 goto low_memory;
2335         natoms = arena->mol->natoms;
2336         size = sizeof(**spp) + sizeof(Vector) * (natoms - 1) * 3;
2337         if (*spp == NULL)
2338                 *spp = (MDSnapshot *)malloc(size);
2339         else
2340                 *spp = (MDSnapshot *)realloc(*spp, size);
2341         if (*spp == NULL)
2342                 goto low_memory;
2343         memset(*spp, 0, size);
2344         (*spp)->step = arena->step;
2345         (*spp)->natoms = natoms;
2346         vp = (*spp)->rvf;
2347         for (i = 0, ap = arena->mol->atoms; i < natoms; i++, ap++) {
2348                 *vp++ = ap->r;
2349                 *vp++ = ap->v;
2350                 *vp++ = ap->f;
2351         }
2352         return;
2353   low_memory:
2354         md_panic(arena, "Low memory while saving snapshot");
2355 }
2356
2357 void
2358 md_restore(MDArena *arena, int idx)
2359 {
2360         int i, natoms, natoms1;
2361         Atom *ap;
2362         Vector *vp;
2363 /*      MDSnapshot **spp; */
2364
2365         if (idx < 0 || idx >= arena->nsnapshots)
2366                 return;
2367         vp = arena->snapshots[idx]->rvf;
2368         natoms = arena->mol->natoms;
2369         natoms1 = arena->snapshots[idx]->natoms;
2370         if (natoms != natoms1) {
2371                 md_warning(arena, "restore: the number of atoms in snapshot (%d) is different from the current number of atoms (%d)\n", natoms1, natoms);
2372         }
2373         for (i = 0, ap = arena->mol->atoms; i < natoms && i < natoms1; i++, ap++) {
2374                 ap->r = *vp++;
2375                 ap->v = *vp++;
2376                 ap->f = *vp++;
2377         }
2378         arena->last_verlet_step = -1;  /*  The Verlet list needs update  */
2379 }
2380
2381 int
2382 md_step(MDArena *arena)
2383 {
2384         md_update_velocities(arena);
2385         md_update_positions(arena);
2386         /* md_rattle_coordinate(arena); */
2387         md_amend_by_symmetry(arena);
2388         calc_force(arena);
2389 /*      md_calc_kinetic_energy(arena); */
2390         md_update_velocities(arena);
2391         /* md_rattle_velocity(arena); */
2392         return 0;
2393 }
2394
2395 void
2396 md_minimize_init(MDArena *arena)
2397 {
2398         int i;
2399         static const Vector zerov = {0, 0, 0};
2400         if (arena->old_forces != NULL)
2401                 free(arena->old_forces);
2402         arena->old_forces = (Vector *)calloc(sizeof(Vector), arena->mol->natoms * 2);
2403         if (arena->old_forces == NULL)
2404                 md_panic(arena, ERROR_out_of_memory);
2405         arena->old_pos = arena->old_forces + arena->mol->natoms;
2406         arena->f_len2 = arena->old_f_len2 = arena->max_gradient = 0.0;
2407         for (i = 0; i < arena->mol->natoms; i++) {
2408                 arena->mol->atoms[i].v = arena->mol->atoms[i].f = zerov;
2409         }
2410         arena->conv_flag = 0;
2411 #if MINIMIZE_CELL
2412         memset(arena->cell_forces, 0, sizeof(Double) * 12);
2413         memset(arena->cell_vels, 0, sizeof(Double) * 12);
2414         memset(arena->old_cell_forces, 0, sizeof(Double) * 12);
2415         memset(arena->old_cell_pars, 0, sizeof(Double) * 12);
2416         arena->cf_len2 = arena->old_cf_len2 = arena->cv_len2 = arena->cell_max_gradient = 0.0;
2417 #endif
2418 }
2419
2420 int
2421 md_minimize_atoms_step(MDArena *arena)
2422 {
2423         Double bk, w1, w2, w3, dump;
2424         Double low, mid, high, low_energy, mid_energy, high_energy, lambda;
2425         Double low_limit, high_limit;
2426         Int i, j, retval, natoms_movable;
2427         Atom *atoms = arena->mol->atoms;
2428         Atom *ap;
2429         Int natoms = arena->mol->natoms;
2430         Vector r, *vp, *vdr;
2431         const Double phi = 0.618033988749895;  /*  The golden ratio  */
2432
2433         md_amend_by_symmetry(arena);
2434
2435         w1 = w2 = 0.0;
2436         retval = 0;
2437         natoms_movable = 0;
2438         for (i = 0, ap = atoms, vp = arena->old_forces; i < natoms; i++, ap++, vp++) {
2439                 if (ap->fix_force < 0)
2440                         continue;
2441                 w1 += VecLength2(ap->f);
2442                 w2 += VecDot(ap->f, *vp);
2443                 natoms_movable++;
2444         }
2445
2446         arena->f_len2 = w1;
2447         if (arena->old_f_len2 == 0.0) {
2448                 /*  New direction  */
2449                 bk = 0.0;
2450         } else {
2451                 bk = (w1 - w2) / arena->old_f_len2;
2452                 if (bk < 0.0)
2453                         bk = 0.0;  /*  New direction  */
2454         }
2455         /*  Update the search direction  */
2456         arena->old_f_len2 = arena->f_len2;
2457         w2 = w3 = 0.0;
2458         dump = 1.0;
2459         for (i = 0, ap = atoms; i < natoms; i++, ap++) {
2460                 if (ap->fix_force < 0)
2461                         continue;
2462                 w1 = VecLength2(ap->f) * dump;
2463                 if (!isfinite(w1))  /*  Is isfinite() available in other platforms?? */
2464                         md_panic(arena, "the gradient at atom %d exceeded limit at step %d", i+1, arena->step);
2465                 if (w1 > 1e4)
2466                         dump *= 1e4 / w1;
2467         }
2468         dump = sqrt(dump);
2469         for (i = 0, ap = atoms, vp = arena->old_forces; i < natoms; i++, ap++, vp++) {
2470                 if (ap->fix_force < 0)
2471                         continue;
2472                 *vp = ap->f;
2473                 *(vp + natoms) = ap->r;
2474                 ap->v.x = ap->v.x * bk + ap->f.x * dump;
2475                 ap->v.y = ap->v.y * bk + ap->f.y * dump;
2476                 ap->v.z = ap->v.z * bk + ap->f.z * dump;
2477                 w1 = VecLength2(ap->v);
2478                 w2 += w1;
2479         /*      if (w1 > 1e4 || !isfinite(w1))
2480                         md_panic(arena, "the gradient at atom %d exceeded limit at step %d", i+1, arena->step); */
2481                 if (w1 > w3)
2482                         w3 = w1;
2483         }
2484         w3 = sqrt(w3);
2485         arena->max_gradient = w3;
2486         arena->v_len2 = w2;
2487 /*      printf("f_len2 = %g, v_len2 = %g, bk = %g, max_grad = %g\n", arena->f_len2, arena->v_len2, bk, w3); */
2488         if (bk == 0.0 && w3 < arena->gradient_convergence)
2489                 return 1;  /*  Gradient is sufficiently small  */
2490
2491         /*  Proceed along ap->v until the energy increases  */
2492         low_limit = arena->coordinate_convergence / arena->max_gradient;
2493         high_limit = 0.1 / arena->max_gradient;
2494         low = 0.0;
2495         low_energy = arena->total_energy;
2496 /*      lambda = 1e-3 * arena->f_len2 / w2; */
2497         lambda = high_limit;
2498         high = lambda;
2499         while (1) {
2500                 for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
2501                         if (ap->fix_force < 0)
2502                                 continue;
2503                         r = ap->r;
2504                         ap->r.x = vp->x + ap->v.x * lambda;
2505                         ap->r.y = vp->y + ap->v.y * lambda;
2506                         ap->r.z = vp->z + ap->v.z * lambda;
2507                         VecDec(r, ap->r);
2508                         VecInc(*vdr, r);
2509                 }
2510                 calc_force(arena);
2511                 mid = lambda;
2512                 mid_energy = arena->total_energy;
2513                 if (mid_energy < low_energy) {
2514                         /*  mid is the 'sufficiently large' step to give lower total energy  */
2515                         if (mid == high) {
2516                                 /*  Higher limit: move by this amount  */
2517                                 retval = 0;
2518                                 goto cleanup;
2519                         }
2520                         break;
2521                 }
2522                 high = mid;
2523                 high_energy = mid_energy;
2524                 lambda *= 0.25;
2525                 if (lambda < low_limit) {
2526                         /*  Cannot find point with lower energy than the starting point  */
2527                         /*  Restore the original position  */
2528                         for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
2529                                 if (ap->fix_force < 0)
2530                                         continue;
2531                                 r = ap->r;
2532                                 ap->r = *vp;
2533                                 VecDec(r, ap->r);
2534                                 VecInc(*vdr, r);
2535                         }
2536                         calc_force(arena);
2537                         lambda = 0.0;
2538                         if (bk == 0.0)
2539                                 retval = 2;  /*  Atom movement is sufficiently small  */
2540                         goto cleanup;
2541                 }
2542         }
2543 /*      printf("Line minimization [%g, %g, %g] (energy [%g, %g, %g]) ", low, mid, high, low_energy, mid_energy, high_energy); */
2544         /*  low_energy >= mid_energy < high_energy  */
2545         /*  Binary search for minimum  */
2546         for (i = 0; i < 5; i++) {
2547                 if (high - mid > mid - low) {
2548                         lambda = high - (high - mid) * phi;
2549                         for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
2550                                 if (ap->fix_force < 0)
2551                                         continue;
2552                                 r = ap->r;
2553                                 ap->r.x = vp->x + ap->v.x * lambda;
2554                                 ap->r.y = vp->y + ap->v.y * lambda;
2555                                 ap->r.z = vp->z + ap->v.z * lambda;
2556                                 VecDec(r, ap->r);
2557                                 VecInc(*vdr, r);
2558                         }
2559                         calc_force(arena);
2560                         if (arena->total_energy < mid_energy) {
2561                                 low = mid;
2562                                 low_energy = mid_energy;
2563                                 mid = lambda;
2564                                 mid_energy = arena->total_energy;
2565                         } else {
2566                                 high = lambda;
2567                                 high_energy = arena->total_energy;
2568                         }
2569                 } else {
2570                         lambda = mid - (mid - low) * phi;
2571                         for (j = 0, ap = atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
2572                                 if (ap->fix_force < 0)
2573                                         continue;
2574                                 r = ap->r;
2575                                 ap->r.x = vp->x + ap->v.x * lambda;
2576                                 ap->r.y = vp->y + ap->v.y * lambda;
2577                                 ap->r.z = vp->z + ap->v.z * lambda;
2578                                 VecDec(r, ap->r);
2579                                 VecInc(*vdr, r);
2580                         }
2581                         calc_force(arena);
2582                         if (arena->total_energy < mid_energy) {
2583                                 high = mid;
2584                                 high_energy = mid_energy;
2585                                 mid = lambda;
2586                                 mid_energy = arena->total_energy;
2587                         } else {
2588                                 low = lambda;
2589                                 low_energy = arena->total_energy;
2590                         }
2591                 }
2592                 if ((high - low) * arena->max_gradient < arena->coordinate_convergence) {
2593                 /*      retval = 2;  *//*  Atom movement is sufficiently small  */
2594                         break;
2595                 }
2596         }
2597   cleanup:
2598 /*      printf("Final lambda = %g (%g)\n", lambda, arena->total_energy); */
2599         return retval;
2600 }
2601
2602 #if 0
2603
2604 static inline void
2605 s_md_modify_atoms_and_cell_parameters(MDArena *arena, Double lambda, Double cell_lambda)
2606 {
2607         Int j, natoms;
2608         Atom *ap;
2609         Vector r, *vp, *vdr;
2610         XtalCell *cell = arena->mol->cell;
2611         natoms = arena->mol->natoms;
2612         for (j = 0, ap = arena->mol->atoms, vp = arena->old_pos, vdr = arena->verlets_dr; j < natoms; j++, ap++, vp++, vdr++) {
2613                 if (ap->fix_force < 0)
2614                         continue;
2615                 r = ap->r;
2616                 ap->r.x = vp->x + ap->v.x * lambda;
2617                 ap->r.y = vp->y + ap->v.y * lambda;
2618                 ap->r.z = vp->z + ap->v.z * lambda;
2619                 VecDec(r, ap->r);
2620                 VecInc(*vdr, r);
2621         }
2622         if (cell_lambda >= 0.0) {
2623                 if (arena->periodic_a) {
2624                         cell->axes[0].x = arena->old_cell_pars[0] + arena->cell_vels[0] * cell_lambda;
2625                         cell->axes[0].y = arena->old_cell_pars[1] + arena->cell_vels[1] * cell_lambda;
2626                         cell->axes[0].z = arena->old_cell_pars[2] + arena->cell_vels[2] * cell_lambda;
2627                 }
2628                 if (arena->periodic_b) {
2629                         cell->axes[1].x = arena->old_cell_pars[3] + arena->cell_vels[3] * cell_lambda;
2630                         cell->axes[1].y = arena->old_cell_pars[4] + arena->cell_vels[4] * cell_lambda;
2631                         cell->axes[1].z = arena->old_cell_pars[5] + arena->cell_vels[5] * cell_lambda;
2632                 }
2633                 if (arena->periodic_c) {
2634                         cell->axes[2].x = arena->old_cell_pars[6] + arena->cell_vels[6] * cell_lambda;
2635                         cell->axes[2].y = arena->old_cell_pars[7] + arena->cell_vels[7] * cell_lambda;
2636                         cell->axes[2].z = arena->old_cell_pars[8] + arena->cell_vels[8] * cell_lambda;
2637                 }
2638                 cell->origin.x = arena->old_cell_pars[9] + arena->cell_vels[9] * cell_lambda;
2639                 cell->origin.y = arena->old_cell_pars[10] + arena->cell_vels[10] * cell_lambda;
2640                 cell->origin.z = arena->old_cell_pars[11] + arena->cell_vels[11] * cell_lambda;
2641                 md_update_cell(arena);
2642         }
2643         md_amend_by_symmetry(arena);
2644 }
2645
2646 int
2647 md_minimize_atoms_and_cell_step(MDArena *arena)
2648 {
2649         Double bk, cbk, w1, w2, w3, w4, damp;
2650         Double low, mid, high, low_energy, mid_energy, high_energy, lambda;
2651         Double low_limit, scale_atom, scale_cell;
2652         Int i, retval, natoms_movable, do_cell;
2653         Atom *atoms = arena->mol->atoms;
2654         Atom *ap;
2655         Int natoms = arena->mol->natoms;
2656         XtalCell *cell;
2657         Vector *vp;
2658         const Double phi = 0.618033988749895;  /*  The golden ratio  */
2659
2660         if (arena->minimize_cell && arena->mol->cell != NULL) {
2661                 do_cell = 1;
2662                 cell = arena->mol->cell;
2663         } else {
2664                 do_cell = 0;
2665                 cell = NULL;
2666         }
2667         bk = cbk = 0.0;
2668
2669         md_amend_by_symmetry(arena);
2670
2671         /*  Get weights of atomic forces  */
2672         w1 = w2 = 0.0;
2673         retval = 0;
2674         natoms_movable = 0;
2675         for (i = 0, ap = atoms, vp = arena->old_forces; i < natoms; i++, ap++, vp++) {
2676                 if (ap->fix_force < 0)
2677                         continue;
2678                 w1 += VecLength2(ap->f);
2679                 w2 += VecDot(ap->f, *vp);
2680                 natoms_movable++;
2681         }
2682         arena->f_len2 = w1;
2683         if (arena->old_f_len2 == 0.0) {
2684                 /*  New direction  */
2685                 bk = 0.0;
2686         } else {
2687                 bk = (w1 - w2) / arena->old_f_len2;
2688                 if (bk < 0.0)
2689                         bk = 0.0;  /*  New direction  */
2690         }
2691
2692         /*  Update the search direction (atoms)  */
2693         arena->old_f_len2 = arena->f_len2;
2694         w2 = w3 = 0.0;
2695         damp = 1.0;
2696         for (i = 0, ap = atoms; i < natoms; i++, ap++) {
2697                 if (ap->fix_force < 0)
2698                         continue;
2699                 w1 = VecLength2(ap->f) * damp;
2700                 if (!isfinite(w1))  /*  Is isfinite() available in other platforms?? */
2701                         md_panic(arena, "the gradient at atom %d exceeded limit at step %d", i+1, arena->step);
2702                 if (w1 > 1e4)
2703                         damp *= 1e4 / w1;
2704         }
2705         damp = sqrt(damp);
2706         for (i = 0, ap = atoms, vp = arena->old_forces; i < natoms; i++, ap++, vp++) {
2707                 if (ap->fix_force < 0)
2708                         continue;
2709                 *vp = ap->f;
2710                 *(vp + natoms) = ap->r;
2711                 ap->v.x = ap->v.x * bk + ap->f.x * damp;
2712                 ap->v.y = ap->v.y * bk + ap->f.y * damp;
2713                 ap->v.z = ap->v.z * bk + ap->f.z * damp;
2714                 w1 = VecLength2(ap->v);
2715                 w2 += w1;
2716                 /*      if (w1 > 1e4 || !isfinite(w1))
2717                  md_panic(arena, "the gradient at atom %d exceeded limit at step %d", i+1, arena->step); */
2718                 if (w1 > w3)
2719                         w3 = w1;
2720         }
2721         w3 = sqrt(w3);
2722         arena->max_gradient = w3;
2723         arena->v_len2 = w2;
2724
2725         /*  Get weights of cell forces  */
2726         if (do_cell) {
2727                 w1 = w2 = 0.0;
2728                 retval = 0;
2729                 for (i = 0; i < 12; i++) {
2730                         w3 = arena->cell_forces[i];
2731                         w1 += w3 * w3;
2732                         w4 = arena->old_cell_forces[i];
2733                         w2 += w3 * w4;
2734                 }
2735                 arena->cf_len2 = w1;
2736                 if (arena->old_cf_len2 == 0.0) {
2737                         /*  New direction  */
2738                         cbk = 0.0;
2739                 } else {
2740                         cbk = (w1 - w2) / arena->old_cf_len2;
2741                         if (cbk < 0.0)
2742                                 cbk = 0.0;  /*  New direction  */
2743                 }
2744
2745                 /*  Update the search direction (cells)  */
2746                 arena->old_cf_len2 = arena->cf_len2;
2747                 damp = 1.0;
2748                 for (i = 0; i < 12; i++) {
2749                         w1 = arena->cell_forces[i];
2750                         w1 = w1 * w1 * damp;
2751                         if (!isfinite(w1))  /*  Is isfinite() available in other platforms?? */
2752                                 md_panic(arena, "the gradient at cell parameter %d exceeded limit at step %d", i+1, arena->step);
2753                         if (w1 > 1e4)
2754                                 damp *= 1e4 / w1;
2755                 }
2756                 damp = sqrt(damp);
2757                 w2 = w3 = 0.0;
2758                 for (i = 0; i < 12; i++) {
2759                         arena->old_cell_forces[i] = w1 = arena->cell_forces[i];
2760                         arena->old_cell_pars[i] = cell->tr[i];
2761                         arena->cell_vels[i] = arena->cell_vels[i] * cbk + w1 * damp;
2762                         w1 *= w1;
2763                         w2 += w1;
2764                         if (w1 > w3)
2765                                 w3 = w1;
2766                 }
2767                 w3 = sqrt(w3);
2768                 arena->cell_max_gradient = w3;
2769                 arena->cv_len2 = w2;
2770         } else {
2771                 cbk = 0.0;
2772                 arena->cell_max_gradient = 0.0;
2773         }
2774
2775         if (bk == 0.0 && arena->max_gradient < arena->gradient_convergence && cbk == 0.0 && arena->cell_max_gradient < arena->gradient_convergence)
2776                 return 1;    /*  Gradient is sufficiently small  */
2777
2778         /*  Proceed along ap->v/cell_vels until the energy increases  */
2779         w1 = arena->coordinate_convergence / arena->max_gradient;
2780         if (do_cell) {
2781                 w2 = arena->coordinate_convergence / arena->cell_max_gradient;
2782                 low_limit = (w1 < w2 ? w1 : w2);
2783         } else low_limit = w1;
2784         scale_atom = 0.1 / arena->max_gradient;
2785         if (do_cell)
2786                 scale_cell = scale_atom;
2787         //      scale_cell = 0.1 / arena->cell_max_gradient;
2788         else scale_cell = -1.0;
2789         low = 0.0;
2790         low_energy = arena->total_energy;
2791         /*      lambda = 1e-3 * arena->f_len2 / w2; */
2792         lambda = 1.0;
2793         high = lambda;
2794         while (1) {
2795                 s_md_modify_atoms_and_cell_parameters(arena, lambda * scale_atom, lambda * scale_cell);
2796                 calc_force(arena);
2797                 mid = lambda;
2798                 mid_energy = arena->total_energy;
2799                 if (mid_energy < low_energy) {
2800                         /*  mid is the 'sufficiently large' step to give lower total energy  */
2801                         if (mid == high) {
2802                                 /*  Higher limit: move by this amount  */
2803                                 retval = 0;
2804                                 goto cleanup;
2805                         }
2806                         break;
2807                 }
2808                 high = mid;
2809                 high_energy = mid_energy;
2810                 lambda *= 0.25;
2811                 if (lambda < low_limit) {
2812                         /*  Cannot find point with lower energy than the starting point  */
2813                         /*  Restore the original position  */
2814                         s_md_modify_atoms_and_cell_parameters(arena, 0, (do_cell ? 0 : -1));
2815                         calc_force(arena);
2816                         lambda = 0.0;
2817                         if (bk == 0.0 && cbk == 0.0)
2818                                 retval = 2;  /*  Movement is sufficiently small  */
2819                         goto cleanup;
2820                 }
2821         }
2822         /*  Binary search for minimum  */
2823         for (i = 0; i < 5; i++) {
2824                 if (high - mid > mid - low) {
2825                         lambda = high - (high - mid) * phi;
2826                         s_md_modify_atoms_and_cell_parameters(arena, lambda * scale_atom, lambda * scale_cell);
2827                         calc_force(arena);
2828                         if (arena->total_energy < mid_energy) {
2829                                 low = mid;
2830                                 low_energy = mid_energy;
2831                                 mid = lambda;
2832                                 mid_energy = arena->total_energy;
2833                         } else {
2834                                 high = lambda;
2835                                 high_energy = arena->total_energy;
2836                         }
2837                 } else {
2838                         lambda = mid - (mid - low) * phi;
2839                         s_md_modify_atoms_and_cell_parameters(arena, lambda * scale_atom, lambda * scale_cell);
2840                         calc_force(arena);
2841                         if (arena->total_energy < mid_energy) {
2842                                 high = mid;
2843                                 high_energy = mid_energy;
2844                                 mid = lambda;
2845                                 mid_energy = arena->total_energy;
2846                         } else {
2847                                 low = lambda;
2848                                 low_energy = arena->total_energy;
2849                         }
2850                 }
2851                 if ((high - low) < low_limit)
2852                         break;
2853         }
2854 cleanup:
2855         /*      printf("Final lambda = %g (%g)\n", lambda, arena->total_energy); */
2856         return retval;
2857 }
2858 #endif
2859
2860
2861 #if MINIMIZE_CELL
2862 static inline void
2863 s_md_modify_cell_parameters(MDArena *arena, Double lambda)
2864 {
2865         XtalCell *cell = arena->mol->cell;
2866         if (arena->periodic_a) {
2867                 cell->axes[0].x = arena->old_cell_pars[0] + arena->cell_vels[0] * lambda;
2868                 cell->axes[0].y = arena->old_cell_pars[1] + arena->cell_vels[1] * lambda;
2869                 cell->axes[0].z = arena->old_cell_pars[2] + arena->cell_vels[2] * lambda;
2870         }
2871         if (arena->periodic_b) {
2872                 cell->axes[1].x = arena->old_cell_pars[3] + arena->cell_vels[3] * lambda;
2873                 cell->axes[1].y = arena->old_cell_pars[4] + arena->cell_vels[4] * lambda;
2874                 cell->axes[1].z = arena->old_cell_pars[5] + arena->cell_vels[5] * lambda;
2875         }
2876         if (arena->periodic_b) {
2877                 cell->axes[2].x = arena->old_cell_pars[6] + arena->cell_vels[6] * lambda;
2878                 cell->axes[2].y = arena->old_cell_pars[7] + arena->cell_vels[7] * lambda;
2879                 cell->axes[2].z = arena->old_cell_pars[8] + arena->cell_vels[8] * lambda;
2880         }
2881         cell->origin.x = arena->old_cell_pars[9] + arena->cell_vels[9] * lambda;
2882         cell->origin.y = arena->old_cell_pars[10] + arena->cell_vels[10] * lambda;
2883         cell->origin.z = arena->old_cell_pars[11] + arena->cell_vels[11] * lambda;
2884         md_update_cell(arena);
2885         md_amend_by_symmetry(arena);
2886 }
2887
2888 int
2889 md_minimize_cell_step(MDArena *arena)
2890 {
2891         Double bk, w1, w2, w3, w4, damp;
2892         Double low, mid, high, low_energy, mid_energy, high_energy, lambda;
2893         Double low_limit, high_limit;
2894         Int i, retval;
2895         XtalCell *cell;
2896         const Double phi = 0.618033988749895;  /*  The golden ratio  */
2897
2898         if (arena->minimize_cell == 0 || (cell = arena->mol->cell) == NULL)
2899                 return 0;
2900
2901         w1 = w2 = 0.0;
2902         retval = 0;
2903         for (i = 0; i < 12; i++) {
2904                 w3 = arena->cell_forces[i];
2905                 w1 += w3 * w3;
2906                 w4 = arena->old_cell_forces[i];
2907                 w2 += w3 * w4;
2908         }
2909
2910         arena->cf_len2 = w1;
2911         if (arena->old_cf_len2 == 0.0) {
2912                 /*  New direction  */
2913                 bk = 0.0;
2914         } else {
2915                 bk = (w1 - w2) / arena->old_cf_len2;
2916                 if (bk < 0.0)
2917                         bk = 0.0;  /*  New direction  */
2918         }
2919         /*  Update the search direction  */
2920         arena->old_cf_len2 = arena->cf_len2;
2921         damp = 1.0;
2922         for (i = 0; i < 12; i++) {
2923                 w1 = arena->cell_forces[i];
2924                 w1 = w1 * w1 * damp;
2925                 if (!isfinite(w1))  /*  Is isfinite() available in other platforms?? */
2926                         md_panic(arena, "the gradient at cell parameter %d exceeded limit at step %d", i+1, arena->step);
2927                 if (w1 > 1e4)
2928                         damp *= 1e4 / w1;
2929         }
2930         damp = sqrt(damp);
2931         w2 = w3 = 0.0;
2932         for (i = 0; i < 12; i++) {
2933                 arena->old_cell_forces[i] = w1 = arena->cell_forces[i];
2934                 arena->old_cell_pars[i] = cell->tr[i];
2935                 arena->cell_vels[i] = arena->cell_vels[i] * bk + w1 * damp;
2936                 w1 *= w1;
2937                 w2 += w1;
2938                 if (w1 > w3)
2939                         w3 = w1;
2940         }
2941         w3 = sqrt(w3);
2942         arena->cell_max_gradient = w3;
2943         arena->cv_len2 = w2;
2944         if (bk == 0.0 && w3 < arena->gradient_convergence)
2945                 return 1;  /*  Gradient is sufficiently small  */
2946
2947         /*  Proceed along cell_vels[] until the energy increases  */
2948         low_limit = arena->coordinate_convergence / arena->max_gradient;
2949         high_limit = 0.1 / arena->max_gradient;
2950         low = 0.0;
2951         low_energy = arena->total_energy;
2952         lambda = high_limit;
2953         high = lambda;
2954         while (1) {
2955                 s_md_modify_cell_parameters(arena, lambda);
2956                 calc_force(arena);
2957                 mid = lambda;
2958                 mid_energy = arena->total_energy;
2959                 if (mid_energy < low_energy) {
2960                         /*  mid is the 'sufficiently large' step to give lower total energy  */
2961                         if (mid == high) {
2962                                 /*  Higher limit: move by this amount  */
2963                                 retval = 0;
2964                                 goto cleanup;
2965                         }
2966                         break;
2967                 }
2968                 high = mid;
2969                 high_energy = mid_energy;
2970                 lambda *= 0.25;
2971                 if (lambda < low_limit) {
2972                         /*  Cannot find point with lower energy than the starting point  */
2973                         /*  Restore the original position  */
2974                         s_md_modify_cell_parameters(arena, 0);
2975                         calc_force(arena);
2976                         lambda = 0.0;
2977                         if (bk == 0.0)
2978                                 retval = 2;  /*  Atom movement is sufficiently small  */
2979                         goto cleanup;
2980                 }
2981         }
2982
2983         /*  low_energy >= mid_energy < high_energy  */
2984         /*  Binary search for minimum  */
2985         for (i = 0; i < 5; i++) {
2986                 if (high - mid > mid - low) {
2987                         lambda = high - (high - mid) * phi;
2988                         s_md_modify_cell_parameters(arena, lambda);
2989                         calc_force(arena);
2990                         if (arena->total_energy < mid_energy) {
2991                                 low = mid;
2992                                 low_energy = mid_energy;
2993                                 mid = lambda;
2994                                 mid_energy = arena->total_energy;
2995                         } else {
2996                                 high = lambda;
2997                                 high_energy = arena->total_energy;
2998                         }
2999                 } else {
3000                         lambda = mid - (mid - low) * phi;
3001                         s_md_modify_cell_parameters(arena, lambda);
3002                         calc_force(arena);
3003                         if (arena->total_energy < mid_energy) {
3004                                 high = mid;
3005                                 high_energy = mid_energy;
3006                                 mid = lambda;
3007                                 mid_energy = arena->total_energy;
3008                         } else {
3009                                 low = lambda;
3010                                 low_energy = arena->total_energy;
3011                         }
3012                 }
3013                 if ((high - low) * arena->cell_max_gradient < arena->coordinate_convergence) {
3014                 /*      retval = 2; */ /*  Atom movement is sufficiently small  */
3015                         break;
3016                 }
3017         }
3018 cleanup:
3019 /*      printf("Cell minimize: ");
3020         for (i = 0; i < 12; i++) {
3021                 printf(" %.6g", arena->cell_vels[i] * lambda);
3022         }
3023         printf("\n"); */
3024
3025         return retval;
3026 }
3027 #endif
3028
3029 int
3030 md_minimize_step(MDArena *arena)
3031 {
3032         int n1, n2;
3033         n1 = md_minimize_atoms_step(arena);
3034         n2 = 0;
3035 #if MINIMIZE_CELL
3036         if (arena->minimize_cell && arena->mol->cell != NULL)
3037                 n2 = md_minimize_cell_step(arena);
3038 #endif
3039         if (n1 > 0 && n2 > 0)
3040                 return n1;
3041         return 0;
3042 }
3043
3044 void
3045 md_update_cell(MDArena *arena)
3046 {
3047         Molecule *mol = arena->mol;
3048         XtalCell *cell = mol->cell;
3049         int i;
3050         MoleculeCalculateCellFromAxes(mol->cell, 1);
3051         for (i = 0; i < mol->nsyms; i++) {
3052                 Transform temp;
3053                 TransformMul(temp, mol->syms[i], cell->rtr);
3054                 TransformMul(arena->cellsyms[i], cell->tr, temp);
3055         }
3056 }
3057
3058 void
3059 md_set_cell(MDArena *arena)
3060 {
3061         Molecule *mol = arena->mol;
3062         if (mol == NULL)
3063                 return;
3064         if (mol->cell != NULL) {
3065                 arena->periodic_a = (mol->cell->flags[0] != 0);
3066                 arena->periodic_b = (mol->cell->flags[1] != 0);
3067                 arena->periodic_c = (mol->cell->flags[2] != 0);
3068                 if (mol->nsyms > 0) {
3069                         Int i;
3070                         if (mol->nsyms != arena->ncellsyms) {
3071                                 arena->ncellsyms = mol->nsyms;
3072                                 arena->cellsyms = (Transform *)realloc(arena->cellsyms, sizeof(Transform) * mol->nsyms);
3073                         }
3074                         for (i = 0; i < mol->nsyms; i++) {
3075                                 Transform temp;
3076                                 TransformMul(temp, mol->syms[i], mol->cell->rtr);
3077                                 TransformMul(arena->cellsyms[i], mol->cell->tr, temp);
3078                         }
3079                 } else {
3080                         arena->ncellsyms = 0;
3081                         if (arena->cellsyms != NULL)
3082                                 free(arena->cellsyms);
3083                         arena->cellsyms = NULL;
3084                 }
3085         } else {
3086                 arena->periodic_a = arena->periodic_b = arena->periodic_c = 0;
3087         }
3088 }
3089
3090 /*  scale_atoms = 1: symmetry-unique atoms are transformed to keep the fractional coordinate constant */
3091 /*  scale_atoms = 2: same as scale_atoms = 1, except that the center of mass of each fragment of connected
3092         atoms are scaled and the atoms within one fragment are moved by the same amount */
3093 /*  NOTE: only symmetry-unique atoms are moved, so don't forget to do md_symmetry_amend()!  */
3094 void
3095 md_scale_cell(MDArena *arena, const Transform tf, int scale_atoms)
3096 {
3097         Transform grad, grad_inv, grad_tr, grad_inv_tr;
3098         int i;
3099         Double volume;
3100         Vector v;
3101         Atom *ap;
3102         XtalCell *cell = arena->mol->cell;
3103         memmove(grad, cell->rtr, sizeof(grad));
3104         VecCross(v, cell->axes[0], cell->axes[1]);
3105         volume = VecDot(v, cell->axes[2]);
3106         TransformVec(&cell->axes[0], tf, &cell->axes[0]);
3107         TransformVec(&cell->axes[1], tf, &cell->axes[1]);
3108         TransformVec(&cell->axes[2], tf, &cell->axes[2]);
3109         md_update_cell(arena);
3110
3111         /*  Deformation gradient (temp) = celltr * old_rcelltr  */
3112         TransformMul(grad, cell->tr, grad);
3113 /*      grad[9] = grad[10] = grad[11] = 0.0;  */
3114         TransformInvert(grad_inv, grad);
3115         memmove(grad_tr, grad, sizeof(grad));
3116         MatrixTranspose(grad_tr, grad_tr);
3117         memmove(grad_inv_tr, grad_inv, sizeof(grad_inv));
3118         MatrixTranspose(grad_inv_tr, grad_inv_tr);
3119
3120         if (scale_atoms == 1) {
3121                 /*  Scale atom positions and velocities  */
3122                 for (i = 0, ap = arena->mol->atoms; i < arena->natoms_uniq; i++, ap++) {
3123                         if (ap->periodic_exclude)
3124                                 continue;
3125                         TransformVec(&ap->r, grad, &ap->r);
3126                         MatrixVec(&ap->f, grad, &ap->f);         /* No translational component */
3127                         MatrixVec(&ap->v, grad_inv_tr, &ap->v);  /* No translational component */
3128                 }
3129         } else if (scale_atoms == 2) {
3130                 int j;
3131                 /*  Scale atom positions and velocities by fragments  */
3132                 for (i = 0; i < arena->nfragments; i++) {
3133                         VecZero(arena->fragment_info[i].pos);
3134                         arena->fragment_info[i].mass = 0.0;
3135                 }
3136                 /*  Get the center of mass for each fragment  */
3137                 for (i = 0, ap = arena->mol->atoms; i < arena->natoms_uniq; i++, ap++) {
3138                         j = arena->fragment_indices[i];
3139                         if (j < 0 || j >= arena->nfragments)
3140                                 continue;
3141                         VecInc(arena->fragment_info[j].pos, ap->r);
3142                         arena->fragment_info[j].mass += ap->weight;
3143                 }
3144                 /*  Calculate the offset for each fragment  */
3145                 for (j = 0; j < arena->nfragments; j++) {
3146                         Vector *vp = &(arena->fragment_info[j].pos);
3147                         if (arena->fragment_info[j].mass > 0.0) {
3148                                 VecScaleSelf(*vp, 1.0 / arena->fragment_info[j].mass);
3149                                 TransformVec(&v, grad, vp);
3150                                 VecSub(*vp, v, *vp);
3151                         } else VecZero(*vp);
3152                 }
3153                 /*  Transform  */
3154                 for (i = 0, ap = arena->mol->atoms; i < arena->natoms_uniq; i++, ap++) {
3155                         if (ap->periodic_exclude)
3156                                 continue;
3157                         MatrixVec(&ap->f, grad, &ap->f);  /*  This is not exact  */
3158                         MatrixVec(&ap->v, grad_inv_tr, &ap->v);
3159                         j = arena->fragment_indices[i];
3160                         if (j < 0 || j >= arena->nfragments)
3161                                 continue;
3162                         VecInc(ap->r, arena->fragment_info[j].pos);
3163                 }
3164         }
3165
3166         arena->last_verlet_step = -1;
3167 }
3168
3169 void
3170 md_flush_output_files(MDArena *arena)
3171 {
3172         if (arena->coord_result != NULL) {
3173                 fflush(arena->coord_result);
3174         }
3175         if (arena->vel_result != NULL) {
3176                 fflush(arena->vel_result);
3177         }
3178         if (arena->force_result != NULL) {
3179                 fflush(arena->force_result);
3180         }
3181         if (arena->extend_result != NULL) {
3182                 fflush(arena->extend_result);
3183         }
3184
3185         if (arena->debug_result != NULL) {
3186                 fflush(arena->debug_result);
3187         }
3188 }
3189
3190 void
3191 md_close_output_files(MDArena *arena)
3192 {
3193         if (arena->coord_result != NULL) {
3194                 fclose(arena->coord_result);
3195                 arena->coord_result = NULL;
3196         }
3197         if (arena->vel_result != NULL) {
3198                 fclose(arena->vel_result);
3199                 arena->vel_result = NULL;
3200         }
3201         if (arena->force_result != NULL) {
3202                 fclose(arena->force_result);
3203                 arena->force_result = NULL;
3204         }
3205         if (arena->extend_result != NULL) {
3206                 fclose(arena->extend_result);
3207                 arena->extend_result = NULL;
3208         }
3209
3210         if (arena->debug_result != NULL) {
3211                 fclose(arena->debug_result);
3212                 arena->debug_result = NULL;
3213         }
3214 }
3215
3216 int
3217 md_copy_coordinates_from_internal(MDArena *arena)
3218 {
3219         /*  Copy the internal r/v/f to xmol  */
3220         int i;
3221         Atom *ap1, *ap2;
3222         if (arena->mol == NULL || arena->xmol == NULL)
3223                 return -1;  /*  Not initialized  */
3224         if (arena->mol->natoms != arena->xmol->natoms)
3225                 return -2;  /*  Number of atoms does not match  */
3226         for (i = 0, ap1 = arena->mol->atoms, ap2 = arena->xmol->atoms; i < arena->mol->natoms; i++, ap1 = ATOM_NEXT(ap1), ap2 = ATOM_NEXT(ap2)) {
3227                 ap2->r = ap1->r;
3228                 ap2->v = ap1->v;
3229                 ap2->f = ap1->f;
3230         }
3231         if (arena->mol->cell != NULL && arena->xmol->cell != NULL)
3232                 memmove(arena->xmol->cell, arena->mol->cell, sizeof(XtalCell));
3233         return 0;
3234 }
3235
3236 int
3237 md_copy_coordinates_to_internal(MDArena *arena)
3238 {
3239         /*  Copy the xmol coordinates to internal  */
3240         int i;
3241         Atom *ap1, *ap2;
3242         if (arena->mol == NULL || arena->xmol == NULL)
3243                 return -1;  /*  Not initialized  */
3244         if (arena->mol->natoms != arena->xmol->natoms)
3245                 return -2;  /*  Number of atoms does not match  */
3246         for (i = 0, ap1 = arena->mol->atoms, ap2 = arena->xmol->atoms; i < arena->mol->natoms; i++, ap1 = ATOM_NEXT(ap1), ap2 = ATOM_NEXT(ap2)) {
3247                 ap1->r = ap2->r;
3248         /*      ap1->occupancy = ap2->occupancy;  *//*  Occupancy can be used to exclude particular atoms  */
3249                 ap1->mm_exclude = ap2->mm_exclude;
3250         }
3251         if (arena->mol->cell != NULL && arena->xmol->cell != NULL)
3252                 memmove(arena->mol->cell, arena->xmol->cell, sizeof(XtalCell));
3253         arena->xmol->needsMDCopyCoordinates = 0;
3254         return 0;
3255 }
3256
3257 int
3258 md_is_running(MDArena *arena)
3259 {
3260         if (arena != NULL && arena->is_running)
3261                 return 1;
3262         else return 0;
3263 }
3264
3265 int
3266 md_main(MDArena *arena, int minimize)
3267 {
3268 //      extern int do_callback(MDArena *);
3269         jmp_buf env;
3270         const char *msg;
3271         int retval = 0;
3272         int (*md_step_func)(MDArena *);
3273
3274         if (arena->is_initialized < 2 || arena->xmol->needsMDRebuild) {
3275                 /*  Prepare MD parameters and runtime fields  */
3276                 msg = md_prepare(arena, 0);
3277                 if (msg != NULL) {
3278                         snprintf(arena->errmsg, sizeof(arena->errmsg), "%s", msg);
3279                         return 1;
3280                 }
3281                 arena->xmol->needsMDCopyCoordinates = 1;  /*  Coordinates will be copied below  */
3282         }
3283
3284         if (arena->xmol->needsMDCopyCoordinates) {
3285                 MoleculeLock(arena->xmol);
3286                 retval = md_copy_coordinates_to_internal(arena);
3287                 MoleculeUnlock(arena->xmol);
3288                 if (retval != 0)
3289                         return retval;
3290                 arena->last_verlet_step = -1;  /*  The Verlet list needs update  */
3291         }
3292
3293         arena->is_running = 1;
3294         arena->is_minimizing = minimize;
3295
3296         if (setjmp(env) == 0) {
3297                 arena->setjmp_buf = &env;
3298                 if (minimize) {
3299                         md_minimize_init(arena);
3300                         md_step_func = md_minimize_step;
3301                         arena->minimize_complete = 0;
3302                 } else {
3303                         md_step_func = md_step;
3304                 }
3305
3306                 /*  Calculate initial energies and forces  */
3307                 arena->step = arena->start_step;
3308                 md_amend_by_symmetry(arena);
3309                 calc_force(arena);
3310                 md_calc_kinetic_energy(arena);
3311                 if (arena->step == 0 || arena->start_step == arena->end_step) {
3312                         md_output_results(arena);
3313                         md_output_energies(arena);
3314                 }
3315                 if (arena->step == 0 && arena->md_callback_func != NULL) {
3316                         retval = (*(arena->md_callback_func))(arena);
3317                         if (retval != 0) {
3318                                 snprintf(arena->errmsg, sizeof(arena->errmsg), "MD Interrupt");
3319                                 goto cleanup;
3320                         }
3321                 }
3322
3323                 /*  Run simulation  */
3324                 for (arena->step = arena->start_step + 1; arena->step <= arena->end_step; arena->step++) {
3325
3326                         /*  Molecules may be modified from the callback procedure  */
3327                         if (arena->is_initialized < 2 || arena->xmol->needsMDRebuild) {
3328                                 msg = md_prepare(arena, 0);
3329                                 if (msg != NULL) {
3330                                         snprintf(arena->errmsg, sizeof(arena->errmsg), "%s", msg);
3331                                         retval = 1;
3332                                         goto cleanup;
3333                                 }
3334                         }
3335
3336                         retval = (*md_step_func)(arena);
3337                         md_calc_kinetic_energy(arena);
3338                         if (!minimize) {
3339                                 if (arena->rescale_temp_freq > 0 && arena->step % arena->rescale_temp_freq == 0)
3340                                         md_scale_velocities(arena);
3341                                 if (arena->reinit_temp_freq > 0 && arena->step % arena->reinit_temp_freq == 0)
3342                                         md_init_velocities(arena);
3343                                 if (arena->andersen_thermo_freq > 0 && arena->step % arena->andersen_thermo_freq == 0)
3344                                         md_andersen_thermostat(arena);
3345                                 if (arena->pressure != NULL)
3346                                         pressure_control(arena);
3347                         }
3348
3349                         if (arena->coord_output_freq > 0 && arena->step % arena->coord_output_freq == 0)
3350                                 md_output_results(arena);
3351                         if (arena->energy_output_freq > 0 && arena->step % arena->energy_output_freq == 0)
3352                                 md_output_energies(arena);
3353
3354                         if (retval != 0) {
3355                                 if (minimize) {
3356                                         switch (retval) {
3357                                                 case 1:
3358                                                         md_log(arena, "Minimize: minimization converged because the maximum gradient becomes less than the threshold.\n");
3359                                                         retval = 0;
3360                                                         arena->minimize_complete = 1;
3361                                                         break;
3362                                                 case 2:
3363                                                         md_log(arena, "Minimize: minimization converged because the maximum movement of atoms becomes less than the threshold.\n");
3364                                                         retval = 0;
3365                                                         arena->minimize_complete = 1;
3366                                                         break;
3367                                         }
3368                                 }
3369                                 goto cleanup;
3370                         }
3371                         if (arena->request_abort != 0) {
3372                                 snprintf(arena->errmsg, sizeof(arena->errmsg), "MD Abort by request");
3373                                 retval = 1;
3374                                 goto cleanup;
3375                         }
3376
3377                         if (arena->md_callback_func != NULL && (arena->callback_freq == 0 || arena->step % arena->callback_freq == 0)) {
3378                                 retval = (*(arena->md_callback_func))(arena);
3379                                 if (retval != 0) {
3380                                         snprintf(arena->errmsg, sizeof(arena->errmsg), "MD Interrupt");
3381                                         goto cleanup;
3382                                 }
3383                         }
3384
3385                 }
3386                 arena->step--;
3387
3388         } else {
3389                 /*  Return from md_panic()  */
3390                 retval = -1;  /*  Some fatal error  */
3391         }
3392
3393 cleanup:
3394         if (retval == 0) {
3395                 arena->start_step = arena->step;  /*  Prepare for next run  */
3396         /*      MoleculeLock(arena->xmol);
3397                 retval = md_copy_coordinates_from_internal(arena);
3398                 MoleculeUnlock(arena->xmol); */
3399         }
3400
3401         arena->setjmp_buf = NULL;
3402
3403         if (arena->is_running) {
3404                 arena->is_running = 0;
3405                 md_flush_output_files(arena);
3406         }
3407
3408         return retval;
3409 }
3410
3411 void
3412 md_set_default(MDArena *arena)
3413 {
3414         arena->start_step = 0;
3415         arena->end_step = 1000;
3416         arena->timestep = 1;
3417
3418         arena->coord_output_freq = 10;
3419         arena->energy_output_freq = 10;
3420         arena->cutoff = 9.0;
3421         arena->electro_cutoff = 9.0;
3422         arena->pairlist_distance = 10.0;
3423         arena->use_xplor_shift = 1;
3424         arena->scale14_vdw = 0.5;
3425         arena->scale14_elect = 0.83;
3426         arena->temperature = 300.0;
3427         arena->velocities_read = 0;
3428         arena->dielectric = 4.8;
3429 /*      arena->probe_radius = 3.2; */
3430         arena->probe_radius = 0.0;
3431         arena->surface_tension = -0.005;
3432         arena->surface_potential_freq = 5;
3433
3434         arena->velocity_limit = 100.0;
3435         arena->sym_tolerance = 5e-3;
3436         arena->gradient_convergence = 1e-6;
3437         arena->coordinate_convergence = 1e-8;
3438         arena->relocate_center = 1;
3439         arena->andersen_thermo_freq = 50;
3440         arena->andersen_thermo_coupling = 0.1;
3441         arena->pressure_freq = 0;
3442
3443         arena->alchem_dlambda = 0.1;
3444
3445 /*      arena->pressure_coupling = 0.4;
3446         arena->pressure_trial_width = 0.01;
3447         arena->pressure_control_algorithm = 0;
3448         arena->pressure_fluctuate_cell_origin = 0.01;
3449         arena->pressure_fluctuate_cell_orientation = 0.01;
3450         arena->pressure[0] = 1;
3451         arena->pressure[1] = 0;
3452         arena->pressure[2] = 0;
3453         arena->pressure[3] = 0;
3454         arena->pressure[4] = 1;
3455         arena->pressure[5] = 0;
3456         arena->pressure[6] = 0;
3457         arena->pressure[7] = 0;
3458         arena->pressure[8] = 1;
3459         arena->cell_flexibility[0] = -1;
3460         arena->cell_flexibility[1] = -1;
3461         arena->cell_flexibility[2] = -1;
3462         arena->cell_flexibility[3] = 0;
3463         arena->cell_flexibility[4] = 0;
3464         arena->cell_flexibility[5] = 0;
3465         arena->cell_flexibility[6] = 0;
3466         arena->cell_flexibility[7] = 0; */
3467 /*      arena->cella.x = 1;
3468         arena->cella.y = 0;
3469         arena->cella.z = 0;
3470         arena->cellb.x = 0;
3471         arena->cellb.y = 1;
3472         arena->cellb.z = 0;
3473         arena->cellc.x = 0;
3474         arena->cellc.y = 0;
3475         arena->cellc.z = 1; */
3476 /*      if (arena->mol != NULL && arena->mol->box != NULL)
3477                 md_update_cell(arena); */
3478 }
3479
3480 MDArena *
3481 md_arena_new(Molecule *xmol)
3482 {
3483         MDArena *arena;
3484         arena = (MDArena *)calloc(sizeof(MDArena), 1);
3485         if (arena == NULL)
3486                 return NULL;
3487         arena->refCount = 1;
3488         md_set_default(arena);
3489         if (xmol != NULL) {
3490                 md_arena_set_molecule(arena, xmol);
3491         }
3492         return arena;
3493 }
3494
3495 MDArena *
3496 md_arena_set_molecule(MDArena *arena, Molecule *xmol)
3497 {
3498         /*  xmol is set to mol  */
3499         if (arena == NULL)
3500                 return NULL;
3501         if (arena->xmol != xmol) {
3502                 if (arena->xmol != NULL) {
3503                         if (arena->xmol->arena == arena)
3504                                 arena->xmol->arena = NULL;
3505                         MoleculeRelease(arena->xmol);
3506                 }
3507                 arena->xmol = xmol;
3508                 if (xmol != NULL) {
3509                         MoleculeRetain(arena->xmol);
3510                         arena->xmol->arena = arena;
3511                 }
3512         }
3513
3514         /*  Dispose the internal cache  */
3515         if (arena->mol != NULL) {
3516                 if (arena->mol->arena == arena)
3517                         arena->mol->arena = NULL;
3518                 MoleculeRelease(arena->mol);
3519                 arena->mol = NULL;
3520         }
3521
3522         if (xmol != NULL) {
3523                 /*  Create an internal copy  */
3524                 Molecule *mol = MoleculeNew();
3525                 Atom *ap, *ap2, arec;
3526                 int i;
3527                 NewArray(&mol->atoms, &mol->natoms, gSizeOfAtomRecord, xmol->natoms);
3528                 for (i = 0, ap = mol->atoms, ap2 = xmol->atoms; i < xmol->natoms; i++, ap = ATOM_NEXT(ap), ap2 = ATOM_NEXT(ap2)) {
3529                         memmove(&arec, ap2, gSizeOfAtomRecord);
3530                         /*  Aniso and frames are unnecessary  */
3531                         arec.aniso = NULL;
3532                         arec.frames = NULL;
3533                         arec.nframes = 0;
3534                         AtomDuplicate(ap, &arec);
3535                 }
3536                 NewArray(&mol->bonds, &mol->nbonds, sizeof(Int) * 2, xmol->nbonds);
3537                 memmove(mol->bonds, xmol->bonds, sizeof(Int) * 2 * xmol->nbonds);
3538                 NewArray(&mol->angles, &mol->nangles, sizeof(Int) * 3, xmol->nangles);
3539                 memmove(mol->angles, xmol->angles, sizeof(Int) * 3 * xmol->nangles);
3540                 NewArray(&mol->dihedrals, &mol->ndihedrals, sizeof(Int) * 4, xmol->ndihedrals);
3541                 memmove(mol->dihedrals, xmol->dihedrals, sizeof(Int) * 4 * xmol->ndihedrals);
3542                 NewArray(&mol->impropers, &mol->nimpropers, sizeof(Int) * 4, xmol->nimpropers);
3543                 memmove(mol->impropers, xmol->impropers, sizeof(Int) * 4 * xmol->nimpropers);
3544                 NewArray(&mol->syms, &mol->nsyms, sizeof(Transform), xmol->nsyms);
3545                 memmove(mol->syms, xmol->syms, sizeof(Transform) * xmol->nsyms);
3546                 if (xmol->cell != NULL) {
3547                         mol->cell = (XtalCell *)malloc(sizeof(XtalCell));
3548                         memmove(mol->cell, xmol->cell, sizeof(XtalCell));
3549                 }
3550                 if (xmol->path != NULL)
3551                         mol->path = strdup(xmol->path);
3552 /*              if (xmol->box != NULL) {
3553                         mol->box = (PeriodicBox *)malloc(sizeof(PeriodicBox));
3554                         memmove(mol->box, xmol->box, sizeof(PeriodicBox));
3555                 } */
3556                 mol->arena = arena;
3557                 mol->par = xmol->par;
3558                 if (mol->par != NULL)
3559                         ParameterRetain(mol->par);
3560                 arena->mol = mol;
3561         }
3562         return arena;
3563 }
3564
3565 MDArena *
3566 md_arena_retain(MDArena *arena)
3567 {
3568         if (arena != NULL)
3569                 arena->refCount++;
3570         return arena;
3571 }
3572
3573 void
3574 md_arena_release(MDArena *arena)
3575 {
3576         int i;
3577         if (arena == NULL)
3578                 return;
3579         if (--arena->refCount != 0)
3580                 return;
3581         if (arena->mol != NULL) {
3582                 if (arena->mol->arena == arena)
3583                         arena->mol->arena = NULL;
3584                 MoleculeRelease(arena->mol);
3585         }
3586         if (arena->xmol != NULL) {
3587                 if (arena->xmol->arena == arena)
3588                         arena->xmol->arena = NULL;
3589                 MoleculeRelease(arena->xmol);
3590         }
3591         if (arena->par != NULL)
3592                 ParameterRelease(arena->par);
3593         if (arena->log_result_name != NULL)
3594                 free((void *)arena->log_result_name);
3595         if (arena->coord_result_name != NULL)
3596                 free((void *)arena->coord_result_name);
3597         if (arena->vel_result_name != NULL)
3598                 free((void *)arena->vel_result_name);
3599         if (arena->force_result_name != NULL)
3600                 free((void *)arena->force_result_name);
3601         if (arena->extend_result_name != NULL)
3602                 free((void *)arena->extend_result_name);
3603         if (arena->debug_result_name != NULL)
3604                 free((void *)arena->debug_result_name);
3605         if (arena->log_result != NULL)
3606                 fclose(arena->log_result);
3607         if (arena->coord_result != NULL)
3608                 fclose(arena->coord_result);
3609         if (arena->vel_result != NULL)
3610                 fclose(arena->vel_result);
3611         if (arena->force_result != NULL)
3612                 fclose(arena->force_result);
3613         if (arena->extend_result != NULL)
3614                 fclose(arena->extend_result);
3615         if (arena->debug_result != NULL)
3616                 fclose(arena->debug_result);
3617         if (arena->custom_bond_pars != NULL)
3618                 free(arena->custom_bond_pars);
3619         if (arena->custom_pars != NULL)
3620                 free(arena->custom_pars);
3621         if (arena->alchem_flags != NULL)
3622                 free(arena->alchem_flags);
3623         if (arena->exforces != NULL)
3624                 free(arena->exforces);
3625         if (arena->bond_par_i != NULL)
3626                 free(arena->bond_par_i);
3627         if (arena->angle_par_i != NULL)
3628                 free(arena->angle_par_i);
3629         if (arena->dihedral_par_i != NULL)
3630                 free(arena->dihedral_par_i);
3631         if (arena->improper_par_i != NULL)
3632                 free(arena->improper_par_i);
3633         if (arena->vdw_par_i != NULL)
3634                 free(arena->vdw_par_i);
3635         if (arena->vdw_cache != NULL)
3636                 free(arena->vdw_cache);
3637         if (arena->energies != NULL)
3638                 free(arena->energies);
3639         if (arena->forces != NULL)
3640                 free(arena->forces);
3641         if (arena->pair_forces != NULL)
3642                 free(arena->pair_forces);
3643         if (arena->cellsyms != NULL)
3644                 free(arena->cellsyms);
3645         if (arena->pressure != NULL)
3646                 pressure_release(arena->pressure);
3647         if (arena->fragment_indices != NULL)
3648                 free(arena->fragment_indices);
3649         if (arena->fragment_info != NULL)
3650                 free(arena->fragment_info);
3651 #warning "TODO: Is arena->exatoms really necessary? "
3652         if (arena->exatoms != NULL)
3653                 free(arena->exatoms);
3654         if (arena->special_positions != NULL)
3655                 free(arena->special_positions);
3656         if (arena->exlist != NULL)
3657                 free(arena->exlist);
3658         if (arena->exinfo != NULL)
3659                 free(arena->exinfo);
3660         if (arena->verlets != NULL)
3661                 free(arena->verlets);
3662         if (arena->verlets_dr != NULL)
3663                 free(arena->verlets_dr);
3664         if (arena->verlet_i != NULL)
3665                 free(arena->verlet_i);
3666         if (arena->snapshots != NULL) {
3667                 for (i = 0; i < arena->nsnapshots; i++) {
3668                         if (arena->snapshots[i] != NULL)
3669                                 free(arena->snapshots[i]);
3670                 }
3671                 free(arena->snapshots);
3672         }
3673         if (arena->old_forces != NULL)
3674                 free(arena->old_forces);
3675         if (arena->old_pos != NULL)
3676                 free(arena->old_pos);
3677         if (arena->graphite != NULL)
3678                 graphite_release(arena->graphite);
3679         if (arena->ring != NULL)
3680                 free(arena->ring);
3681         free(arena);
3682 }
3683
3684 void
3685 md_arena_init_from_arena(MDArena *arena, MDArena *another_arena)
3686 {
3687         if (arena == NULL || another_arena == NULL)
3688                 return;
3689
3690         arena->timestep = another_arena->timestep;
3691
3692         arena->coord_output_freq = another_arena->coord_output_freq;
3693         arena->energy_output_freq = another_arena->energy_output_freq;
3694         arena->cutoff = another_arena->cutoff;
3695         arena->electro_cutoff = another_arena->electro_cutoff;
3696         arena->pairlist_distance = another_arena->pairlist_distance;
3697         arena->use_xplor_shift = another_arena->use_xplor_shift;
3698         arena->scale14_vdw = another_arena->scale14_vdw;
3699         arena->scale14_elect = another_arena->scale14_elect;
3700         arena->temperature = another_arena->temperature;
3701         arena->rescale_temp_freq = another_arena->rescale_temp_freq;
3702         arena->reinit_temp_freq = another_arena->reinit_temp_freq;
3703         arena->velocities_read = another_arena->velocities_read;
3704         arena->dielectric = another_arena->dielectric;
3705         arena->probe_radius = another_arena->probe_radius;
3706         arena->surface_tension = another_arena->surface_tension;
3707         arena->surface_potential_freq = another_arena->surface_potential_freq;
3708         arena->anbond_thres = another_arena->anbond_thres;
3709         arena->anbond_anneal_rate = another_arena->anbond_anneal_rate;
3710         arena->velocity_limit = another_arena->velocity_limit;
3711         arena->sym_tolerance = another_arena->sym_tolerance;
3712         arena->gradient_convergence = another_arena->gradient_convergence;
3713         arena->coordinate_convergence = another_arena->coordinate_convergence;
3714         arena->relocate_center = another_arena->relocate_center;
3715         arena->andersen_thermo_freq = another_arena->andersen_thermo_freq;
3716         arena->andersen_thermo_coupling = another_arena->andersen_thermo_coupling;
3717         arena->pressure_freq = another_arena->pressure_freq;
3718 }
3719
3720 /*
3721 PROGRAM
3722 {
3723         call force(f)
3724         do loop=1,nstep
3725           v = v + (dt/2) * (f/m)
3726           r = r + dt*v;
3727           call rattle_coodinate
3728           call force(f)
3729           v = v + (dt/2) * (f/m)
3730           call rattle_velocity
3731     enddo
3732 }
3733 */