blob: d866b4d3350f621bdbb8a09e8a0fdb7b60ffe3c4 [file] [log] [blame]
/*************************************************************************/
/* */
/* Copyright (c) 1994 Stanford University */
/* */
/* All rights reserved. */
/* */
/* Permission is given to use, copy, and modify this software for any */
/* non-commercial purpose as long as this copyright notice is not */
/* removed. All other uses, including redistribution in whole or in */
/* part, are forbidden without prior written permission. */
/* */
/* This software is provided with absolutely no warranty and no */
/* support. */
/* */
/*************************************************************************/
EXTERN_ENV
#include "mdvar.h"
#include "parameters.h"
#include "mddata.h"
#include "split.h"
#include "global.h"
/* predicts new values for displacement and its five derivatives
*
* NOR1 : NOR1 = NORDER + 1 = 7 (for a sixth-order method)
*/
void PREDIC(double *C, long NOR1, long ProcID)
{
/* this routine calculates predicted F(X), F'(X), F''(X), ... */
long JIZ;
long JI;
long L;
double S;
long func, mol, dir, atom;
JIZ=2;
/* .....loop over F(X), F'(X), F''(X), ..... */
for (func = 0; func < NORDER; func++) {
for (mol = StartMol[ProcID]; mol < StartMol[ProcID+1]; mol++)
for ( dir = 0; dir < NDIR; dir++)
for ( atom = 0; atom < NATOM; atom++ ) {
JI = JIZ;
/* sum over Taylor Series */
S = 0.0;
for ( L = func; L < NORDER; L++) {
S += C[JI] * VAR[mol].F[L+1][dir][atom];
JI++;
} /* for */
VAR[mol].F[func][dir][atom] += S;
} /* for atom */
JIZ += NOR1;
} /* for func */
} /* end of subroutine PREDIC */
/* corrects the predicted values, based on forces etc. computed in the interim
*
* PCC : the predictor-corrector constants
* NOR1: NORDER + 1 = 7 for a sixth-order method)
*/
void CORREC(double *PCC, long NOR1, long ProcID)
{
/*
.....this routine calculates corrected F(X), F'(X), F"(X), ....
from corrected F(X) = predicted F(X) + PCC(1)*(FR-SD)
where SD is predicted accl. F"(X) and FR is computed
accl. (force/mass) at predicted position
*/
double Y;
long mol, dir, atom, func;
for (mol = StartMol[ProcID]; mol < StartMol[ProcID+1]; mol++) {
for (dir = 0; dir < NDIR; dir++) {
for (atom = 0; atom < NATOM; atom++) {
Y = VAR[mol].F[FORCES][dir][atom] - VAR[mol].F[ACC][dir][atom];
for ( func = 0; func < NOR1; func++)
VAR[mol].F[func][dir][atom] += PCC[func] * Y;
} /* for atom */
} /* for dir */
} /* for mol */
} /* end of subroutine CORREC */