/* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "stdlib.h" #include "string.h" #include "stdio.h" #include "fix_rattle.h" #include "atom.h" #include "atom_vec.h" #include "molecule.h" #include "update.h" #include "respa.h" #include "modify.h" #include "domain.h" #include "force.h" #include "bond.h" #include "angle.h" #include "comm.h" #include "group.h" #include "fix_respa.h" #include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; /* ---------------------------------------------------------------------- */ FixRattle::FixRattle(LAMMPS *lmp, int narg, char **arg) : FixShake(lmp, narg, arg) { // define timesteps with conversion factors dt = update->dt; dtfvo2 = update->dt * force->ftm2v/2.; dtfsq_full = update->dt * update->dt * force->ftm2v; // allocate memory for unconstrained velocity update vp = NULL; grow_arrays(atom->nmax); // DEBUG: process id and stream for debugging purposes MPI_Comm_rank(MPI_COMM_WORLD, &pid); dfp = (RATTLE_DEBUG) ? stdout : NULL; } FixRattle::~FixRattle() { memory->destroy(vp); } /**************************************************************************** * * Purpose: This method carries out the unconstrained velocity update * and iterates over all clusters (molecules). * * NOTE: This function modifies the velocities v[i][j] and vp[i][j]. * *****************************************************************************/ void FixRattle::post_force(int vflag) { // DEBUG if (RATTLE_DEBUG) fprintf(dfp,"FixRattle(%d)::post_force: head of function.\n", update->ntimestep); // this method carries out an unconstrained velocity update by half a timestep (similar to FixShake::unconstrained_update()) update_v_half_nocons(); // communicate the unconstrained velocities if (nprocs > 1) { comm_mode = VP; comm->forward_comm_fix(this); } // correct the velocity for each molecule accordingly int m; for (int i = 0; i < nlist; i++) { m = list[i]; if (shake_flag[m] == 2) vrattle(m); // not implemented yet else if (shake_flag[m] == 3) vrattle3(m); // not implemented yet else if (shake_flag[m] == 4) vrattle4(m); // not implemented yet else { vrattle3angle(m); } } // DEBUG if (RATTLE_DEBUG) fprintf(dfp, "Proc(%d), FixRattle(%d)::post_force: At the end of the function.\n", pid, update->ntimestep); } /**************************************************************************** * * Purpose: This function removes the velocity components along the bonds * by using constraining forces such that the velocities * satisfy v_ij(t+dt) * r_ij(t+dt) = 0 for i,j in molecule m. * * NOTE: This function modifies the velocities v[i][j] and vp[i][j]. * The above constraints will only be satisfied after * fix_nve::final_integrate() was executed. * *****************************************************************************/ void FixRattle::vrattle3angle(int m) { // 1.) setup coordinate and velocity vectors for molecule m int nlist,list[3]; double imass[3]; double c[3], l[3], a[3][3]; // local atom IDs and constraint distances tagint i0 = atom->map(shake_atom[m][0]); tagint i1 = atom->map(shake_atom[m][1]); tagint i2 = atom->map(shake_atom[m][2]); double bond1 = bond_distance[shake_type[m][0]]; double bond2 = bond_distance[shake_type[m][1]]; double bond12 = angle_distance[shake_type[m][2]]; // r01,r02,r12 = distance vec between atoms, with PBC double r01[3]; r01[0] = x[i1][0] - x[i0][0]; r01[1] = x[i1][1] - x[i0][1]; r01[2] = x[i1][2] - x[i0][2]; domain->minimum_image(r01); double r02[3]; r02[0] = x[i2][0] - x[i0][0]; r02[1] = x[i2][1] - x[i0][1]; r02[2] = x[i2][2] - x[i0][2]; domain->minimum_image(r02); double r12[3]; r12[0] = x[i2][0] - x[i1][0]; r12[1] = x[i2][1] - x[i1][1]; r12[2] = x[i2][2] - x[i1][2]; domain->minimum_image(r12); // v01, v02, v12 = distance vecs for velocities after unconstrained update double vp01[3]; vp01[0] = vp[i1][0] - vp[i0][0]; vp01[1] = vp[i1][1] - vp[i0][1]; vp01[2] = vp[i1][2] - vp[i0][2]; double vp02[3]; vp02[0] = vp[i2][0] - vp[i0][0]; vp02[1] = vp[i2][1] - vp[i0][1]; vp02[2] = vp[i2][2] - vp[i0][2]; double vp12[3]; vp12[0] = vp[i2][0] - vp[i1][0]; vp12[1] = vp[i2][1] - vp[i1][1]; vp12[2] = vp[i2][2] - vp[i1][2]; // scalar distances between atoms double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; double r12sq = r12[0]*r12[0] + r12[1]*r12[1] + r12[2]*r12[2]; double r01dotr02 = r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]; double r01dotr12 = r01[0]*r12[0] + r01[1]*r12[1] + r01[2]*r12[2]; double r02dotr12 = r02[0]*r12[0] + r02[1]*r12[1] + r02[2]*r12[2]; // matrix coeffs and rhs for lamda equations if (rmass) { imass[0] = 1.0/rmass[i0]; imass[1] = 1.0/rmass[i1]; imass[2] = 1.0/rmass[i2]; } else { imass[0] = 1.0/mass[type[i0]]; imass[1] = 1.0/mass[type[i1]]; imass[2] = 1.0/mass[type[i2]]; } //2.) Calculate Lagrange parameteres // setup matrix a[0][0] = dt * (imass[1] + imass[0]) * r01sq; a[0][1] = dt * (imass[0] ) * r01dotr02; a[0][2] = dt * (-imass[1] ) * r01dotr12; a[1][0] = a[0][1]; a[1][1] = dt * (imass[0] + imass[2]) * r02sq; a[1][2] = dt * (imass[2] ) * r02dotr12; a[2][0] = a[0][2]; a[2][1] = a[1][2]; a[2][2] = dt * (imass[2] + imass[1]) * r12sq; // setup RHS c[0] = - (vp01[0]*r01[0] + vp01[1]*r01[1] + vp01[2]*r01[2]); c[1] = - (vp02[0]*r02[0] + vp02[1]*r02[1] + vp02[2]*r02[2]); c[2] = - (vp12[0]*r12[0] + vp12[1]*r12[1] + vp12[2]*r12[2]); // calculate the inverse matrix with a method of choice, a := inverse(a) solve3x3exactly(a,c,l); //3.) add corrections to the velocities if the process owns this atom if (i0 < nlocal) { for (int k=0; k<3; k++) v[i0][k] -= dt * imass[0] * ( l[0] * r01[k] + l[1] * r02[k] ) ; } if (i1 < nlocal) { for (int k=0; k<3; k++) v[i1][k] -= dt * imass[1] * ( -l[0] * r01[k] + l[2] * r12[k] ) ; } if (i2 < nlocal) { for (int k=0; k<3; k++) v[i2][k] -= dt * imass[2] * ( -l[1] * r02[k] - l[2] * r12[k] ) ; } } void FixRattle::solve3x3exactly(const double a[][3], const double c[], double l[]) { double ai[3][3]; double determ, determinv; // calculate the determinant of the matrix determ = a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + a[0][2]*a[1][0]*a[2][1] - a[0][0]*a[1][2]*a[2][1] - a[0][1]*a[1][0]*a[2][2] - a[0][2]*a[1][1]*a[2][0]; // check if matrix is actually invertible if (determ == 0.0) error->one(FLERR,"Rattle determinant = 0.0"); // calculate the inverse 3x3 matrix: A^(-1) = (ai_jk) determinv = 1.0/determ; ai[0][0] = determinv * (a[1][1]*a[2][2] - a[1][2]*a[2][1]); ai[0][1] = -determinv * (a[0][1]*a[2][2] - a[0][2]*a[2][1]); ai[0][2] = determinv * (a[0][1]*a[1][2] - a[0][2]*a[1][1]); ai[1][0] = -determinv * (a[1][0]*a[2][2] - a[1][2]*a[2][0]); ai[1][1] = determinv * (a[0][0]*a[2][2] - a[0][2]*a[2][0]); ai[1][2] = -determinv * (a[0][0]*a[1][2] - a[0][2]*a[1][0]); ai[2][0] = determinv * (a[1][0]*a[2][1] - a[1][1]*a[2][0]); ai[2][1] = -determinv * (a[0][0]*a[2][1] - a[0][1]*a[2][0]); ai[2][2] = determinv * (a[0][0]*a[1][1] - a[0][1]*a[1][0]); // Calcualte the solution: (l01, l02, l12)^T = A^(-1) * c for (int i=0; i<3; i++) { l[i] = 0; for (int j=0; j<3; j++) l[i] += ai[i][j] * c[j]; } } // TODO: Implement void FixRattle::vrattle(int m) { } void FixRattle::vrattle3(int m) { } void FixRattle::vrattle4(int m) { } int FixRattle::setmask() { int mask = 0; mask |= PRE_NEIGHBOR; mask |= POST_FORCE; mask |= INITIAL_INTEGRATE; // DEBUG if (RATTLE_TEST) mask |= END_OF_STEP; return mask; } /**************************************************************************** * * Purpose: Add constraining forces to f[i][j] such that the constraints * on the coordinates will be satisfied after * fix_nve::initial_integrate() was called. * * The constrining forces are calculated such that the conditions * * (r_ij(t+dt) - d_ij)^2 = 0 for i,j in molecule m, * * holds, where d_ij are the bond lengths to be imposed. * * *****************************************************************************/ void FixRattle::initial_integrate(int vflag) { // change timestep for compatibility with FixShake implementation FixShake::dtfsq = dtfsq_full/2.; // QUESTION: Where should this be? if (update->ntimestep == next_output) stats(); // unconstrained coordinate update (xshake) FixShake::unconstrained_update(); // communicate unconstrained coordinates (xshake) if running in parallel if (nprocs > 1) { comm_mode = XSHAKE; comm->forward_comm_fix(this); } // QUESTION: Where should this be? My suggestion is FixRattle::post_force // FixShake::shake3angle updates the virial if evflag = 1. // Should I temporarily set evflag to 0 and then in FixRattle::post_force // back to its initial status? if (vflag) v_setup(vflag); else evflag = 0; // loop over clusters to add constraint forces int m; for (int i = 0; i < nlist; i++) { m = list[i]; if (shake_flag[m] == 2) shake(m); else if (shake_flag[m] == 3) shake3(m); else if (shake_flag[m] == 4) shake4(m); else shake3angle(m); } // restore the timestep again FixShake::dtfsq = dtfsq_full; } /**************************************************************************** * * Purpose: Integrate the velocities by half a timestep and store * result in vp[i][j]. * *****************************************************************************/ void FixRattle::update_v_half_nocons() { // Just for debugging if (RATTLE_DEBUG) fprintf(dfp,"FixRattle(%d)::update_v_half_nocons: unconstrained velocity update.\n", update->ntimestep); double dtfvinvm; if (rmass) { for (int i = 0; i < nlocal; i++) { if (shake_flag[i]) { dtfvinvm = dtfvo2 / rmass[i]; for (int k=0; k<3; k++) vp[i][k] = v[i][k] + dtfvinvm * f[i][k]; } else vp[i][0] = vp[i][1] = vp[i][2] = 0; } } else { for (int i = 0; i < nlocal; i++) { dtfvinvm = dtfvo2/mass[type[i]]; if (shake_flag[i]) { for (int k=0; k<3; k++) vp[i][k] = v[i][k] + dtfvinvm * f[i][k]; } else vp[i][0] = vp[i][1] = vp[i][2] = 0; } } } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double FixRattle::memory_usage() { int nmax = atom->nmax; double bytes = FixShake::memory_usage(); bytes += nmax*3 * sizeof(double); return bytes; } /* ---------------------------------------------------------------------- allocate local atom-based arrays ------------------------------------------------------------------------- */ void FixRattle::grow_arrays(int nmax) { FixShake::grow_arrays(nmax); memory->destroy(vp); memory->create(vp,nmax,3,"rattle:vp"); } int FixRattle::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; double dx,dy,dz; m = 0; switch (comm_mode) { case XSHAKE: // SHAKE packs xshake m = FixShake::pack_forward_comm(n, list, buf, pbc_flag, pbc); break; case VP: for (i = 0; i < n; i++) { j = list[i]; buf[m++] = vp[j][0]; buf[m++] = vp[j][1]; buf[m++] = vp[j][2]; } break; // DEBUG ONLY case V: for (i = 0; i < n; i++) { j = list[i]; buf[m++] = v[j][0]; buf[m++] = v[j][1]; buf[m++] = v[j][2]; } break; default: printf("Proc(%d), FixRattle::pack_forward_comm(%d): ERROR: Ended up in default branch of switch(comm_mode)!\n", pid, update->ntimestep); exit(1); } return m; } void FixRattle::unpack_forward_comm(int n, int first, double *buf) { int i, m, last; m = 0; last = first + n; switch (comm_mode) { case XSHAKE: // SHAKE packs xshake FixShake::unpack_forward_comm(n, first,buf); break; case VP: for (i = first; i < last; i++) { vp[i][0] = buf[m++]; vp[i][1] = buf[m++]; vp[i][2] = buf[m++]; } break; // DEBUG ONLY case V: for (i = first; i < last; i++) { v[i][0] = buf[m++]; v[i][1] = buf[m++]; v[i][2] = buf[m++]; } break; default: printf("Proc(%d), FixRattle::unpack_forward_comm(%d): ERROR: Ended up in default branch of switch(comm_mode)! comm_mode = %d\n", pid, update->ntimestep, comm_mode); exit(1); } } /********************************* * * FOR DEBUGGING ONLY * *********************************/ void FixRattle::end_of_step() { if (RATTLE_DEBUG) fprintf(dfp,"FixRattle(%d)::end_of_step: head of function.\n", update->ntimestep); // communicate xshake and vp if (nprocs > 1) { comm_mode = V; comm->forward_comm_fix(this); } if (!check_constraints(v, true, true)) { printf( "FixRattle(%d)::end_of_step: ERROR: RATTLE DOES NOT WORK CORRECTLY!\n", update->ntimestep); exit(1); } else if (RATTLE_DEBUG) fprintf(dfp, "FixRattle(%d)::end_of_step: coordinates and velocities are OK!\n", update->ntimestep); } bool FixRattle::check_molecule(double **v, int m, bool checkr, bool checkv) { // local atom IDs and constraint distances int i0 = atom->map(shake_atom[m][0]); int i1 = atom->map(shake_atom[m][1]); int i2 = atom->map(shake_atom[m][2]); double bond1 = bond_distance[shake_type[m][0]]; double bond2 = bond_distance[shake_type[m][1]]; double bond12 = angle_distance[shake_type[m][2]]; // r01,r02,r12 = distance vec between atoms, with PBC double r01[3]; r01[0] = x[i1][0] - x[i0][0]; r01[1] = x[i1][1] - x[i0][1]; r01[2] = x[i1][2] - x[i0][2]; domain->minimum_image(r01); double r02[3]; r02[0] = x[i2][0] - x[i0][0]; r02[1] = x[i2][1] - x[i0][1]; r02[2] = x[i2][2] - x[i0][2]; domain->minimum_image(r02); double r12[3]; r12[0] = x[i2][0] - x[i1][0]; r12[1] = x[i2][1] - x[i1][1]; r12[2] = x[i2][2] - x[i1][2]; domain->minimum_image(r12); // v01, v02, v12 = distance vecs for velocities after unconstrained update double v01[3]; v01[0] = v[i1][0] - v[i0][0]; v01[1] = v[i1][1] - v[i0][1]; v01[2] = v[i1][2] - v[i0][2]; double v02[3]; v02[0] = v[i2][0] - v[i0][0]; v02[1] = v[i2][1] - v[i0][1]; v02[2] = v[i2][2] - v[i0][2]; double v12[3]; v12[0] = v[i2][0] - v[i1][0]; v12[1] = v[i2][1] - v[i1][1]; v12[2] = v[i2][2] - v[i1][2]; // scalar distances between atoms double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2]; double r02sq = r02[0]*r02[0] + r02[1]*r02[1] + r02[2]*r02[2]; double r12sq = r12[0]*r12[0] + r12[1]*r12[1] + r12[2]*r12[2]; double r01v01 = r01[0]*v01[0] + r01[1]*v01[1] + r01[2]*v01[2]; double r02v02 = r02[0]*v02[0] + r02[1]*v02[1] + r02[2]*v02[2]; double r12v12 = r12[0]*v12[0] + r12[1]*v12[1] + r12[2]*v12[2]; // check velocity constraints bool stat = true; const double tol = tolerance; if (checkr && ( fabs(sqrt(r01sq) - bond1) > tol || fabs(sqrt(r02sq) - bond2) > tol || fabs(sqrt(r12sq) - bond12) > tol ) ) { printf("\n\n\n"); printf(" ************* POSITION ERROR ***************\n\n"); printf("MOLECULE:\t%d\n", m); printf("LOCAL INDICES:\t%d %d %d\n", i0,i1,i2); printf("GLOBAL INDICES:\t%d %d %d\n", atom->tag[i0], atom->tag[i1], atom->tag[i2]); printf("BOND 1 (0<->1): %lf\n", bond1); printf("BOND 2 (0<->2): %lf\n", bond2); printf("BOND 3 (1<->2): %lf\n", bond12); printf("r01: %lf\n", sqrt(r01sq)); printf("r02: %lf\n", sqrt(r02sq)); printf("r12: %lf\n", sqrt(r12sq)); printf(" --> Coordinate constrains failed!\n"); stat = false; } if (checkv && ( fabs(r01v01) > tol || fabs(r02v02) > tol || fabs(r12v12) > tol) ) { printf("\n\n************* VELOCITY ERROR ***************\n"); printf("MOLECULE:\t%d\n", m); printf("LOCAL INDICES:\t%d %d %d\n", i0,i1,i2); printf("GLOBAL INDICES:\t%d %d %d\n", atom->tag[i0], atom->tag[i1], atom->tag[i2]); printf("< r01,v01 > = %lf\n", r01v01); printf("< r02,v02 > = %lf\n", r02v02); printf("< r12,v12 > = %lf\n", r12v12); printf(" --> Velocity constraints failed!\n"); stat = false; } return stat; } bool FixRattle::check_constraints(double **v, bool checkr, bool checkv) { int m; bool stat = true; int i=0; while (i < nlist && stat) { m = list[i]; stat = stat && check_molecule(v, m, checkr, checkv); i++; } return stat; }