mirror of https://github.com/lammps/lammps.git
2291 lines
70 KiB
C++
2291 lines
70 KiB
C++
/* ----------------------------------------------------------------------
|
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
www.cs.sandia.gov/~sjplimp/lammps.html
|
|
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
|
|
|
|
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_shake.h"
|
|
#include "atom.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 "memory.h"
|
|
#include "error.h"
|
|
|
|
#define BIG 1.0e20
|
|
#define MASSDELTA 0.1
|
|
#define MIN(a,b) ((a) < (b) ? (a) : (b))
|
|
#define MAX(a,b) ((a) > (b) ? (a) : (b))
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
FixShake::FixShake(int narg, char **arg) : Fix(narg, arg)
|
|
{
|
|
MPI_Comm_rank(world,&me);
|
|
MPI_Comm_size(world,&nprocs);
|
|
|
|
PI = 4.0*atan(1.0);
|
|
|
|
// error check
|
|
|
|
if (atom->molecular == 0)
|
|
error->all("Cannot use fix shake with non-molecular system");
|
|
|
|
// perform initial allocation of atom-based arrays
|
|
// register with atom class
|
|
|
|
shake_flag = NULL;
|
|
shake_atom = shake_type = NULL;
|
|
xshake = NULL;
|
|
|
|
grow_arrays(atom->nmax);
|
|
atom->add_callback(0);
|
|
|
|
// parse SHAKE args
|
|
|
|
if (narg < 8) error->all("Illegal fix shake command");
|
|
|
|
tolerance = atof(arg[3]);
|
|
max_iter = atoi(arg[4]);
|
|
output_every = atoi(arg[5]);
|
|
|
|
// parse SHAKE args for bond and angle types
|
|
// will be used by find_clusters
|
|
// store args for "b" "a" "t" as flags in (1:n) list for fast access
|
|
// store args for "m" in list of length nmass for looping over
|
|
// for "m" verify that atom masses have been set
|
|
|
|
bond_flag = new int[atom->nbondtypes+1];
|
|
for (int i = 1; i <= atom->nbondtypes; i++) bond_flag[i] = 0;
|
|
angle_flag = new int[atom->nangletypes+1];
|
|
for (int i = 1; i <= atom->nangletypes; i++) angle_flag[i] = 0;
|
|
type_flag = new int[atom->ntypes+1];
|
|
for (int i = 1; i <= atom->ntypes; i++) type_flag[i] = 0;
|
|
mass_list = new double[atom->ntypes];
|
|
nmass = 0;
|
|
|
|
char mode = '\0';
|
|
int next = 6;
|
|
while (next < narg) {
|
|
|
|
if (strcmp(arg[next],"b") == 0) mode = 'b';
|
|
else if (strcmp(arg[next],"a") == 0) mode = 'a';
|
|
else if (strcmp(arg[next],"t") == 0) mode = 't';
|
|
else if (strcmp(arg[next],"m") == 0) {
|
|
mode = 'm';
|
|
atom->check_mass();
|
|
|
|
} else if (mode == 'b') {
|
|
int i = atoi(arg[next]);
|
|
if (i < 1 || i > atom->nbondtypes)
|
|
error->all("Invalid bond type index for fix shake");
|
|
bond_flag[i] = 1;
|
|
|
|
} else if (mode == 'a') {
|
|
int i = atoi(arg[next]);
|
|
if (i < 1 || i > atom->nangletypes)
|
|
error->all("Invalid angle type index for fix shake");
|
|
angle_flag[i] = 1;
|
|
|
|
} else if (mode == 't') {
|
|
int i = atoi(arg[next]);
|
|
if (i < 1 || i > atom->ntypes)
|
|
error->all("Invalid atom type index for fix shake");
|
|
type_flag[i] = 1;
|
|
|
|
} else if (mode == 'm') {
|
|
double rmass = atof(arg[next]);
|
|
if (rmass == 0.0) error->all("Invalid atom mass for fix shake");
|
|
if (nmass == atom->ntypes) error->all("Too many masses for fix shake");
|
|
mass_list[nmass++] = rmass;
|
|
|
|
} else error->all("Illegal fix shake command");
|
|
next++;
|
|
}
|
|
|
|
// allocate bond and angle distance arrays, indexed from 1 to n
|
|
|
|
bond_distance = new double[atom->nbondtypes+1];
|
|
angle_distance = new double[atom->nangletypes+1];
|
|
|
|
// allocate statistics arrays
|
|
|
|
if (output_every) {
|
|
int nb = atom->nbondtypes + 1;
|
|
b_count = new int[nb];
|
|
b_count_all = new int[nb];
|
|
b_ave = new double[nb];
|
|
b_ave_all = new double[nb];
|
|
b_max = new double[nb];
|
|
b_max_all = new double[nb];
|
|
b_min = new double[nb];
|
|
b_min_all = new double[nb];
|
|
|
|
int na = atom->nangletypes + 1;
|
|
a_count = new int[na];
|
|
a_count_all = new int[na];
|
|
a_ave = new double[na];
|
|
a_ave_all = new double[na];
|
|
a_max = new double[na];
|
|
a_max_all = new double[na];
|
|
a_min = new double[na];
|
|
a_min_all = new double[na];
|
|
}
|
|
|
|
// identify all SHAKE clusters
|
|
|
|
find_clusters();
|
|
|
|
// initialize list of SHAKE clusters to constrain
|
|
|
|
maxlist = 0;
|
|
list = NULL;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
FixShake::~FixShake()
|
|
{
|
|
// if atom class still exists:
|
|
// unregister this fix so atom class doesn't invoke it any more
|
|
// set bond_type and angle_type back to positive for SHAKE clusters
|
|
// must set for all SHAKE bonds and angles stored by each atom
|
|
|
|
if (atom) {
|
|
atom->delete_callback(id,0);
|
|
|
|
int **bond_type = atom->bond_type;
|
|
int **angle_type = atom->angle_type;
|
|
int nlocal = atom->nlocal;
|
|
|
|
int n;
|
|
for (int i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i] == 0) continue;
|
|
else if (shake_flag[i] == 1) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = anglefind(i,shake_atom[i][1],shake_atom[i][2]);
|
|
if (n >= 0) angle_type[i][n] = -angle_type[i][n];
|
|
} else if (shake_flag[i] == 2) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
} else if (shake_flag[i] == 3) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
} else if (shake_flag[i] == 4) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][3]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
}
|
|
}
|
|
}
|
|
|
|
// delete locally stored arrays
|
|
|
|
memory->sfree(shake_flag);
|
|
memory->destroy_2d_int_array(shake_atom);
|
|
memory->destroy_2d_int_array(shake_type);
|
|
memory->destroy_2d_double_array(xshake);
|
|
|
|
delete [] bond_flag;
|
|
delete [] angle_flag;
|
|
delete [] type_flag;
|
|
delete [] mass_list;
|
|
|
|
delete [] bond_distance;
|
|
delete [] angle_distance;
|
|
|
|
if (output_every) {
|
|
delete [] b_count;
|
|
delete [] b_count_all;
|
|
delete [] b_ave;
|
|
delete [] b_ave_all;
|
|
delete [] b_max;
|
|
delete [] b_max_all;
|
|
delete [] b_min;
|
|
delete [] b_min_all;
|
|
|
|
delete [] a_count;
|
|
delete [] a_count_all;
|
|
delete [] a_ave;
|
|
delete [] a_ave_all;
|
|
delete [] a_max;
|
|
delete [] a_max_all;
|
|
delete [] a_min;
|
|
delete [] a_min_all;
|
|
}
|
|
|
|
memory->sfree(list);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
int FixShake::setmask()
|
|
{
|
|
int mask = 0;
|
|
mask |= PRE_NEIGHBOR;
|
|
mask |= POST_FORCE;
|
|
mask |= POST_FORCE_RESPA;
|
|
return mask;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
set bond and angle distances
|
|
this init must happen after force->bond and force->angle inits
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::init()
|
|
{
|
|
int i,m,flag,flag_all,type1,type2,bond1_type,bond2_type;
|
|
double rsq,angle;
|
|
|
|
// error if more than one shake fix
|
|
|
|
int count = 0;
|
|
for (i = 0; i < modify->nfix; i++)
|
|
if (strcmp(modify->fix[i]->style,"shake") == 0) count++;
|
|
if (count > 1) error->all("More than one shake fix");
|
|
|
|
// error if npt,nph fix comes before shake fix
|
|
|
|
for (i = 0; i < modify->nfix; i++) {
|
|
if (strcmp(modify->fix[i]->style,"npt") == 0) break;
|
|
if (strcmp(modify->fix[i]->style,"nph") == 0) break;
|
|
}
|
|
if (i < modify->nfix) {
|
|
for (int j = i; j < modify->nfix; j++)
|
|
if (strcmp(modify->fix[j]->style,"shake") == 0)
|
|
error->all("Shake fix must come before NPT/NPH fix");
|
|
}
|
|
|
|
// if rRESPA, find associated fix that must exist
|
|
// could have changed locations in fix list since created
|
|
// set ptrs to rRESPA variables
|
|
|
|
if (strcmp(update->integrate_style,"respa") == 0) {
|
|
for (i = 0; i < modify->nfix; i++)
|
|
if (strcmp(modify->fix[i]->style,"RESPA") == 0) ifix_respa = i;
|
|
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
|
loop_respa = ((Respa *) update->integrate)->loop;
|
|
step_respa = ((Respa *) update->integrate)->step;
|
|
}
|
|
|
|
// set communication size in comm class
|
|
|
|
comm->maxforward_fix = MAX(comm->maxforward_fix,3);
|
|
|
|
// set equilibrium bond distances
|
|
|
|
if (force->bond == NULL)
|
|
error->all("Bond potential must be defined for SHAKE");
|
|
for (i = 1; i <= atom->nbondtypes; i++)
|
|
bond_distance[i] = force->bond->equilibrium_distance(i);
|
|
|
|
// set equilibrium angle distances
|
|
|
|
int nlocal = atom->nlocal;
|
|
|
|
for (i = 1; i <= atom->nangletypes; i++) {
|
|
if (angle_flag[i] == 0) continue;
|
|
|
|
// scan all atoms for a SHAKE angle cluster
|
|
// extract bond types for the 2 bonds in the cluster
|
|
// bond types must be same in all clusters of this angle type,
|
|
// else set error flag
|
|
|
|
flag = 0;
|
|
bond1_type = bond2_type = 0;
|
|
for (m = 0; m < nlocal; m++) {
|
|
if (shake_flag[m] != 1) continue;
|
|
if (shake_type[m][2] != i) continue;
|
|
type1 = MIN(shake_type[m][0],shake_type[m][1]);
|
|
type2 = MAX(shake_type[m][0],shake_type[m][1]);
|
|
if (bond1_type > 0) {
|
|
if (type1 != bond1_type || type2 != bond2_type) {
|
|
flag = 1;
|
|
break;
|
|
}
|
|
}
|
|
bond1_type = type1;
|
|
bond2_type = type2;
|
|
}
|
|
|
|
// error check for any bond types that are not the same
|
|
|
|
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world);
|
|
if (flag_all) error->all("Shake angles have different bond types");
|
|
|
|
// insure all procs have bond types
|
|
|
|
MPI_Allreduce(&bond1_type,&flag_all,1,MPI_INT,MPI_MAX,world);
|
|
bond1_type = flag_all;
|
|
MPI_Allreduce(&bond2_type,&flag_all,1,MPI_INT,MPI_MAX,world);
|
|
bond2_type = flag_all;
|
|
|
|
// if bond types are 0, no SHAKE angles of this type exist
|
|
// just skip this angle
|
|
|
|
if (bond1_type == 0) {
|
|
angle_distance[i] = 0.0;
|
|
continue;
|
|
}
|
|
|
|
// compute the angle distance as a function of 2 bond distances
|
|
|
|
angle = force->angle->equilibrium_angle(i);
|
|
rsq = 2.0*bond_distance[bond1_type]*bond_distance[bond2_type] *
|
|
(1.0-cos(angle));
|
|
angle_distance[i] = sqrt(rsq);
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
SHAKE as pre-integrator constraint
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::setup()
|
|
{
|
|
pre_neighbor();
|
|
|
|
if (output_every) stats();
|
|
|
|
// setup SHAKE output
|
|
|
|
int ntimestep = update->ntimestep;
|
|
next_output = ntimestep + output_every;
|
|
if (output_every == 0) next_output = update->laststep + 1;
|
|
if (output_every && ntimestep % output_every != 0)
|
|
next_output = (ntimestep/output_every)*output_every + output_every;
|
|
|
|
// half timestep constraint on pre-step, full timestep thereafter
|
|
|
|
if (strcmp(update->integrate_style,"verlet") == 0) {
|
|
dtv = update->dt;
|
|
dtfsq = 0.5 * update->dt * update->dt * force->ftm2v;
|
|
post_force(1);
|
|
dtfsq = update->dt * update->dt * force->ftm2v;
|
|
} else {
|
|
dtv = step_respa[0];
|
|
dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v;
|
|
dtf_inner = dtf_innerhalf;
|
|
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
|
|
post_force_respa(1,nlevels_respa-1,0);
|
|
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
|
|
dtf_inner = step_respa[0] * force->ftm2v;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
build list of SHAKE clusters to constrain
|
|
if one or more atoms in cluster are on this proc,
|
|
this proc lists the cluster exactly once
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::pre_neighbor()
|
|
{
|
|
int atom1,atom2,atom3,atom4;
|
|
|
|
// local copies of atom quantities
|
|
// used by SHAKE until next re-neighboring
|
|
|
|
x = atom->x;
|
|
v = atom->v;
|
|
f = atom->f;
|
|
mass = atom->mass;
|
|
type = atom->type;
|
|
nlocal = atom->nlocal;
|
|
|
|
// extend size of SHAKE list if necessary
|
|
|
|
if (nlocal > maxlist) {
|
|
maxlist = nlocal;
|
|
memory->sfree(list);
|
|
list = (int *) memory->smalloc(maxlist*sizeof(int),"shake:list");
|
|
}
|
|
|
|
// build list of SHAKE clusters I compute
|
|
|
|
nlist = 0;
|
|
|
|
for (int i = 0; i < nlocal; i++)
|
|
if (shake_flag[i]) {
|
|
if (shake_flag[i] == 2) {
|
|
atom1 = atom->map(shake_atom[i][0]);
|
|
atom2 = atom->map(shake_atom[i][1]);
|
|
if (atom1 == -1 || atom2 == -1) {
|
|
char str[128];
|
|
sprintf(str,"Shake atoms %d %d missing on proc %d at step %d",
|
|
shake_atom[i][0],shake_atom[i][1],me,update->ntimestep);
|
|
error->one(str);
|
|
}
|
|
if (i <= atom1 && i <= atom2) list[nlist++] = i;
|
|
} else if (shake_flag[i] % 2 == 1) {
|
|
atom1 = atom->map(shake_atom[i][0]);
|
|
atom2 = atom->map(shake_atom[i][1]);
|
|
atom3 = atom->map(shake_atom[i][2]);
|
|
if (atom1 == -1 || atom2 == -1 || atom3 == -1) {
|
|
char str[128];
|
|
sprintf(str,"Shake atoms %d %d %d missing on proc %d at step %d",
|
|
shake_atom[i][0],shake_atom[i][1],shake_atom[i][2],
|
|
me,update->ntimestep);
|
|
error->one(str);
|
|
}
|
|
if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
|
|
} else {
|
|
atom1 = atom->map(shake_atom[i][0]);
|
|
atom2 = atom->map(shake_atom[i][1]);
|
|
atom3 = atom->map(shake_atom[i][2]);
|
|
atom4 = atom->map(shake_atom[i][3]);
|
|
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) {
|
|
char str[128];
|
|
sprintf(str,"Shake atoms %d %d %d %d missing on proc %d at step %d",
|
|
shake_atom[i][0],shake_atom[i][1],
|
|
shake_atom[i][2],shake_atom[i][3],
|
|
me,update->ntimestep);
|
|
error->one(str);
|
|
}
|
|
if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
|
|
list[nlist++] = i;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
compute the force adjustment for SHAKE constraint
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::post_force(int vflag_in)
|
|
{
|
|
if (update->ntimestep == next_output) stats();
|
|
|
|
// xshake = unconstrained move with current v,f
|
|
|
|
unconstrained_update();
|
|
|
|
// communicate results if necessary
|
|
|
|
if (nprocs > 1) comm->comm_fix(this);
|
|
|
|
// zero out SHAKE contribution to virial
|
|
|
|
vflag = vflag_in;
|
|
if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0;
|
|
|
|
// loop over clusters
|
|
|
|
int m;
|
|
for (int i = 0; i < nlist; i++) {
|
|
m = list[i];
|
|
if (shake_flag[m] == 2) shake2(m);
|
|
else if (shake_flag[m] == 3) shake3(m);
|
|
else if (shake_flag[m] == 4) shake4(m);
|
|
else shake3angle(m);
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
count # of degrees-of-freedom removed by SHAKE for atoms in igroup
|
|
------------------------------------------------------------------------- */
|
|
|
|
int FixShake::dof(int igroup)
|
|
{
|
|
int groupbit = group->bitmask[igroup];
|
|
|
|
int *mask = atom->mask;
|
|
int *tag = atom->tag;
|
|
int nlocal = atom->nlocal;
|
|
|
|
// count dof in a cluster if and only if the central atom is in group
|
|
// and this atom i is the central atom
|
|
|
|
int n = 0;
|
|
for (int i = 0; i < nlocal; i++) {
|
|
if (!(mask[i] & groupbit)) continue;
|
|
if (shake_flag[i] == 0) continue;
|
|
if (shake_atom[i][0] != tag[i]) continue;
|
|
if (shake_flag[i] == 1) n += 3;
|
|
else if (shake_flag[i] == 2) n += 1;
|
|
else if (shake_flag[i] == 3) n += 2;
|
|
else if (shake_flag[i] == 4) n += 3;
|
|
}
|
|
|
|
int nall;
|
|
MPI_Allreduce(&n,&nall,1,MPI_INT,MPI_SUM,world);
|
|
return nall;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
identify whether each atom is in a SHAKE cluster
|
|
only include atoms in fix group and those bonds/angles specified in input
|
|
test whether all clusters are valid
|
|
set shake_flag, shake_atom, shake_type values
|
|
set bond,angle types negative so will be ignored in neighbor lists
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::find_clusters()
|
|
{
|
|
int i,j,m,n;
|
|
int flag,flag_all,messtag,loop,nbuf,nbufmax,size;
|
|
double imass,jmass,rmass;
|
|
int *buf,*bufcopy;
|
|
MPI_Request request;
|
|
MPI_Status status;
|
|
|
|
if (me == 0 && screen) fprintf(screen,"Finding SHAKE clusters ...\n");
|
|
|
|
// local copies of atom ptrs
|
|
|
|
int *tag = atom->tag;
|
|
int *type = atom->type;
|
|
int *mask = atom->mask;
|
|
double *mass = atom->mass;
|
|
int **bond_type = atom->bond_type;
|
|
int **angle_type = atom->angle_type;
|
|
int **nspecial = atom->nspecial;
|
|
int **special = atom->special;
|
|
int nlocal = atom->nlocal;
|
|
|
|
// setup ring of procs
|
|
|
|
int next = me + 1;
|
|
int prev = me -1;
|
|
if (next == nprocs) next = 0;
|
|
if (prev < 0) prev = nprocs - 1;
|
|
|
|
// -----------------------------------------------------
|
|
// allocate arrays for self (1d) and bond partners (2d)
|
|
// max = max # of bond partners for owned atoms = 2nd dim of partner arrays
|
|
// npartner[i] = # of bonds attached to atom i
|
|
// nshake[i] = # of SHAKE bonds attached to atom i
|
|
// partner_tag[i][] = global IDs of each partner
|
|
// partner_mask[i][] = mask of each partner
|
|
// partner_type[i][] = type of each partner
|
|
// partner_bondtype[i][] = type of bond attached to each partner
|
|
// partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not
|
|
// partner_nshake[i][] = nshake value for each partner
|
|
// -----------------------------------------------------
|
|
|
|
int max = 0;
|
|
for (i = 0; i < nlocal; i++) max = MAX(max,nspecial[i][0]);
|
|
|
|
int *npartner = (int *)
|
|
memory->smalloc(nlocal*sizeof(double),"shake:npartner");
|
|
int *nshake = (int *)
|
|
memory->smalloc(nlocal*sizeof(double),"shake:nshake");
|
|
int **partner_tag =
|
|
memory->create_2d_int_array(nlocal,max,"shake:partner_tag");
|
|
int **partner_mask =
|
|
memory->create_2d_int_array(nlocal,max,"shake:partner_mask");
|
|
int **partner_type =
|
|
memory->create_2d_int_array(nlocal,max,"shake:partner_type");
|
|
int **partner_bondtype =
|
|
memory->create_2d_int_array(nlocal,max,"shake:partner_bondtype");
|
|
int **partner_shake =
|
|
memory->create_2d_int_array(nlocal,max,"shake:partner_shake");
|
|
int **partner_nshake =
|
|
memory->create_2d_int_array(nlocal,max,"shake:partner_nshake");
|
|
|
|
// -----------------------------------------------------
|
|
// set npartner and partner_tag from special arrays
|
|
// -----------------------------------------------------
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
npartner[i] = nspecial[i][0];
|
|
for (j = 0; j < npartner[i]; j++) partner_tag[i][j] = special[i][j];
|
|
}
|
|
|
|
// -----------------------------------------------------
|
|
// set partner_mask, partner_type, partner_bondtype for bonded partners
|
|
// requires communication for off-proc partners
|
|
// -----------------------------------------------------
|
|
|
|
// fill in mask, type, bondtype if own bond partner
|
|
// info to store in buf for each off-proc bond =
|
|
// 2 atoms IDs in bond, space for mask, type, bondtype
|
|
// nbufmax = largest buffer needed to hold info from any proc
|
|
|
|
nbuf = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
partner_mask[i][j] = 0;
|
|
partner_type[i][j] = 0;
|
|
partner_bondtype[i][j] = 0;
|
|
|
|
m = atom->map(partner_tag[i][j]);
|
|
if (m >= 0 && m < nlocal) {
|
|
partner_mask[i][j] = mask[m];
|
|
partner_type[i][j] = type[m];
|
|
n = bondfind(i,tag[i],partner_tag[i][j]);
|
|
if (n >= 0) partner_bondtype[i][j] = bond_type[i][n];
|
|
else {
|
|
n = bondfind(m,tag[i],partner_tag[i][j]);
|
|
if (n >= 0) partner_bondtype[i][j] = bond_type[m][n];
|
|
}
|
|
} else nbuf += 5;
|
|
}
|
|
}
|
|
|
|
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
|
|
|
|
buf = new int[nbufmax];
|
|
bufcopy = new int[nbufmax];
|
|
|
|
// fill buffer with info
|
|
|
|
size = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
m = atom->map(partner_tag[i][j]);
|
|
if (m < 0 || m >= nlocal) {
|
|
buf[size] = tag[i];
|
|
buf[size+1] = partner_tag[i][j];
|
|
buf[size+2] = 0;
|
|
buf[size+3] = 0;
|
|
n = bondfind(i,tag[i],partner_tag[i][j]);
|
|
if (n >= 0) buf[size+4] = bond_type[i][n];
|
|
else buf[size+4] = 0;
|
|
size += 5;
|
|
}
|
|
}
|
|
}
|
|
|
|
// cycle buffer around ring of procs back to self
|
|
// when receive buffer, scan bond partner IDs for atoms I own
|
|
// if I own partner, fill in mask and type, search for bond with 1st atom
|
|
// and fill in bondtype
|
|
|
|
messtag = 1;
|
|
for (loop = 0; loop < nprocs; loop++) {
|
|
i = 0;
|
|
while (i < size) {
|
|
m = atom->map(buf[i+1]);
|
|
if (m >= 0 && m < nlocal) {
|
|
buf[i+2] = mask[m];
|
|
buf[i+3] = type[m];
|
|
if (buf[i+4] == 0) {
|
|
n = bondfind(m,buf[i],buf[i+1]);
|
|
if (n >= 0) buf[i+4] = bond_type[m][n];
|
|
}
|
|
}
|
|
i += 5;
|
|
}
|
|
if (me != next) {
|
|
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
|
|
MPI_Send(buf,size,MPI_INT,next,messtag,world);
|
|
MPI_Wait(&request,&status);
|
|
MPI_Get_count(&status,MPI_INT,&size);
|
|
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
|
|
}
|
|
}
|
|
|
|
// store partner info returned to me
|
|
|
|
m = 0;
|
|
while (m < size) {
|
|
i = atom->map(buf[m]);
|
|
for (j = 0; j < npartner[i]; j++)
|
|
if (buf[m+1] == partner_tag[i][j]) break;
|
|
partner_mask[i][j] = buf[m+2];
|
|
partner_type[i][j] = buf[m+3];
|
|
partner_bondtype[i][j] = buf[m+4];
|
|
m += 5;
|
|
}
|
|
|
|
delete [] buf;
|
|
delete [] bufcopy;
|
|
|
|
// error check for unfilled partner info
|
|
// if partner_type not set, is an error
|
|
// partner_bondtype may not be set if special list is not consistent
|
|
// with bondatom (e.g. due to delete_bonds command)
|
|
// this is OK if one or both atoms are not in fix group, since
|
|
// bond won't be SHAKEn anyway
|
|
// else it's an error
|
|
|
|
flag = 0;
|
|
for (i = 0; i < nlocal; i++)
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
if (partner_type[i][j] == 0) flag = 1;
|
|
if (!(mask[i] & groupbit)) continue;
|
|
if (!(partner_mask[i][j] & groupbit)) continue;
|
|
if (partner_bondtype[i][j] == 0) flag = 1;
|
|
}
|
|
|
|
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
|
|
if (flag_all) error->all("Did not find fix shake partner info");
|
|
|
|
// -----------------------------------------------------
|
|
// identify SHAKEable bonds
|
|
// set nshake[i] = # of SHAKE bonds attached to atom i
|
|
// set partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not
|
|
// both atoms must be in group, bondtype must be > 0
|
|
// check if bondtype is in input bond_flag
|
|
// check if type of either atom is in input type_flag
|
|
// check if mass of either atom is in input mass_list
|
|
// -----------------------------------------------------
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
nshake[i] = 0;
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
partner_shake[i][j] = 0;
|
|
|
|
if (!(mask[i] & groupbit)) continue;
|
|
if (!(partner_mask[i][j] & groupbit)) continue;
|
|
if (partner_bondtype[i][j] <= 0) continue;
|
|
|
|
if (bond_flag[partner_bondtype[i][j]]) {
|
|
partner_shake[i][j] = 1;
|
|
nshake[i]++;
|
|
continue;
|
|
}
|
|
if (type_flag[type[i]] || type_flag[partner_type[i][j]]) {
|
|
partner_shake[i][j] = 1;
|
|
nshake[i]++;
|
|
continue;
|
|
}
|
|
if (nmass) {
|
|
imass = mass[type[i]];
|
|
jmass = mass[partner_type[i][j]];
|
|
for (m = 0; m < nmass; m++) {
|
|
rmass = mass_list[m];
|
|
if (fabs(rmass-imass) <= MASSDELTA ||
|
|
fabs(rmass-jmass) <= MASSDELTA) {
|
|
partner_shake[i][j] = 1;
|
|
nshake[i]++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// -----------------------------------------------------
|
|
// set partner_nshake for bonded partners
|
|
// requires communication for off-proc partners
|
|
// -----------------------------------------------------
|
|
|
|
// fill in partner_nshake if own bond partner
|
|
// info to store in buf for each off-proc bond =
|
|
// 2 atoms IDs in bond, space for nshake value
|
|
// nbufmax = largest buffer needed to hold info from any proc
|
|
|
|
nbuf = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
m = atom->map(partner_tag[i][j]);
|
|
if (m >= 0 && m < nlocal) partner_nshake[i][j] = nshake[m];
|
|
else nbuf += 3;
|
|
}
|
|
}
|
|
|
|
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
|
|
|
|
buf = new int[nbufmax];
|
|
bufcopy = new int[nbufmax];
|
|
|
|
// fill buffer with info
|
|
|
|
size = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
m = atom->map(partner_tag[i][j]);
|
|
if (m < 0 || m >= nlocal) {
|
|
buf[size] = tag[i];
|
|
buf[size+1] = partner_tag[i][j];
|
|
size += 3;
|
|
}
|
|
}
|
|
}
|
|
|
|
// cycle buffer around ring of procs back to self
|
|
// when receive buffer, scan bond partner IDs for atoms I own
|
|
// if I own partner, fill in nshake value
|
|
|
|
messtag = 2;
|
|
for (loop = 0; loop < nprocs; loop++) {
|
|
i = 0;
|
|
while (i < size) {
|
|
m = atom->map(buf[i+1]);
|
|
if (m >= 0 && m < nlocal) buf[i+2] = nshake[m];
|
|
i += 3;
|
|
}
|
|
if (me != next) {
|
|
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
|
|
MPI_Send(buf,size,MPI_INT,next,messtag,world);
|
|
MPI_Wait(&request,&status);
|
|
MPI_Get_count(&status,MPI_INT,&size);
|
|
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
|
|
}
|
|
}
|
|
|
|
// store partner info returned to me
|
|
|
|
m = 0;
|
|
while (m < size) {
|
|
i = atom->map(buf[m]);
|
|
for (j = 0; j < npartner[i]; j++)
|
|
if (buf[m+1] == partner_tag[i][j]) break;
|
|
partner_nshake[i][j] = buf[m+2];
|
|
m += 3;
|
|
}
|
|
|
|
delete [] buf;
|
|
delete [] bufcopy;
|
|
|
|
// -----------------------------------------------------
|
|
// error checks
|
|
// no atom with nshake > 3
|
|
// no connected atoms which both have nshake > 1
|
|
// -----------------------------------------------------
|
|
|
|
flag = 0;
|
|
for (i = 0; i < nlocal; i++) if (nshake[i] > 3) flag = 1;
|
|
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
|
|
if (flag_all) error->all("Shake cluster of more than 4 atoms");
|
|
|
|
flag = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
if (nshake[i] <= 1) continue;
|
|
for (j = 0; j < npartner[i]; j++)
|
|
if (partner_shake[i][j] && partner_nshake[i][j] > 1) flag = 1;
|
|
}
|
|
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
|
|
if (flag_all) error->all("Shake clusters are connected");
|
|
|
|
// -----------------------------------------------------
|
|
// set SHAKE arrays that are stored with atoms & add angle constraints
|
|
// zero shake arrays for all owned atoms
|
|
// if I am central atom set shake_flag & shake_atom & shake_type
|
|
// for 2-atom clusters, I am central atom if my atom ID < partner ID
|
|
// for 3-atom clusters, test for angle constraint
|
|
// angle will be stored by this atom if it exists
|
|
// if angle type matches angle_flag, then it is angle-constrained
|
|
// shake_flag[] = 0 if atom not in SHAKE cluster
|
|
// 2,3,4 = size of bond-only cluster
|
|
// 1 = 3-atom angle cluster
|
|
// shake_atom[][] = global IDs of 2,3,4 atoms in cluster
|
|
// central atom is 1st
|
|
// for 2-atom cluster, lowest ID is 1st
|
|
// shake_type[][] = bondtype of each bond in cluster
|
|
// for 3-atom angle cluster, 3rd value is angletype
|
|
// -----------------------------------------------------
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
shake_flag[i] = 0;
|
|
shake_atom[i][0] = 0;
|
|
shake_atom[i][1] = 0;
|
|
shake_atom[i][2] = 0;
|
|
shake_atom[i][3] = 0;
|
|
shake_type[i][0] = 0;
|
|
shake_type[i][1] = 0;
|
|
shake_type[i][2] = 0;
|
|
|
|
if (nshake[i] == 1) {
|
|
for (j = 0; j < npartner[i]; j++)
|
|
if (partner_shake[i][j]) break;
|
|
if (partner_nshake[i][j] == 1 && tag[i] < partner_tag[i][j]) {
|
|
shake_flag[i] = 2;
|
|
shake_atom[i][0] = tag[i];
|
|
shake_atom[i][1] = partner_tag[i][j];
|
|
shake_type[i][0] = partner_bondtype[i][j];
|
|
}
|
|
}
|
|
|
|
if (nshake[i] > 1) {
|
|
shake_flag[i] = 1;
|
|
shake_atom[i][0] = tag[i];
|
|
for (j = 0; j < npartner[i]; j++)
|
|
if (partner_shake[i][j]) {
|
|
m = shake_flag[i];
|
|
shake_atom[i][m] = partner_tag[i][j];
|
|
shake_type[i][m-1] = partner_bondtype[i][j];
|
|
shake_flag[i]++;
|
|
}
|
|
}
|
|
|
|
if (nshake[i] == 2) {
|
|
n = anglefind(i,shake_atom[i][1],shake_atom[i][2]);
|
|
if (n < 0) continue;
|
|
if (angle_type[i][n] < 0) continue;
|
|
if (angle_flag[angle_type[i][n]]) {
|
|
shake_flag[i] = 1;
|
|
shake_type[i][2] = angle_type[i][n];
|
|
}
|
|
}
|
|
}
|
|
|
|
// -----------------------------------------------------
|
|
// set shake_flag,shake_atom,shake_type for non-central atoms
|
|
// requires communication for off-proc atoms
|
|
// -----------------------------------------------------
|
|
|
|
// fill in shake arrays for each bond partner I own
|
|
// info to store in buf for each off-proc bond =
|
|
// all values from shake_flag, shake_atom, shake_type
|
|
// nbufmax = largest buffer needed to hold info from any proc
|
|
|
|
nbuf = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i] == 0) continue;
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
if (partner_shake[i][j] == 0) continue;
|
|
m = atom->map(partner_tag[i][j]);
|
|
if (m >= 0 && m < nlocal) {
|
|
shake_flag[m] = shake_flag[i];
|
|
shake_atom[m][0] = shake_atom[i][0];
|
|
shake_atom[m][1] = shake_atom[i][1];
|
|
shake_atom[m][2] = shake_atom[i][2];
|
|
shake_atom[m][3] = shake_atom[i][3];
|
|
shake_type[m][0] = shake_type[i][0];
|
|
shake_type[m][1] = shake_type[i][1];
|
|
shake_type[m][2] = shake_type[i][2];
|
|
} else nbuf += 9;
|
|
}
|
|
}
|
|
|
|
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
|
|
|
|
buf = new int[nbufmax];
|
|
bufcopy = new int[nbufmax];
|
|
|
|
// fill buffer with info
|
|
|
|
size = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i] == 0) continue;
|
|
for (j = 0; j < npartner[i]; j++) {
|
|
if (partner_shake[i][j] == 0) continue;
|
|
m = atom->map(partner_tag[i][j]);
|
|
if (m < 0 || m >= nlocal) {
|
|
buf[size] = partner_tag[i][j];
|
|
buf[size+1] = shake_flag[i];
|
|
buf[size+2] = shake_atom[i][0];
|
|
buf[size+3] = shake_atom[i][1];
|
|
buf[size+4] = shake_atom[i][2];
|
|
buf[size+5] = shake_atom[i][3];
|
|
buf[size+6] = shake_type[i][0];
|
|
buf[size+7] = shake_type[i][1];
|
|
buf[size+8] = shake_type[i][2];
|
|
size += 9;
|
|
}
|
|
}
|
|
}
|
|
|
|
// cycle buffer around ring of procs back to self
|
|
// when receive buffer, scan for ID that I own
|
|
// if I own ID, fill in shake array values
|
|
|
|
messtag = 3;
|
|
for (loop = 0; loop < nprocs; loop++) {
|
|
i = 0;
|
|
while (i < size) {
|
|
m = atom->map(buf[i]);
|
|
if (m >= 0 && m < nlocal) {
|
|
shake_flag[m] = buf[i+1];
|
|
shake_atom[m][0] = buf[i+2];
|
|
shake_atom[m][1] = buf[i+3];
|
|
shake_atom[m][2] = buf[i+4];
|
|
shake_atom[m][3] = buf[i+5];
|
|
shake_type[m][0] = buf[i+6];
|
|
shake_type[m][1] = buf[i+7];
|
|
shake_type[m][2] = buf[i+8];
|
|
}
|
|
i += 9;
|
|
}
|
|
if (me != next) {
|
|
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
|
|
MPI_Send(buf,size,MPI_INT,next,messtag,world);
|
|
MPI_Wait(&request,&status);
|
|
MPI_Get_count(&status,MPI_INT,&size);
|
|
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
|
|
}
|
|
}
|
|
|
|
delete [] buf;
|
|
delete [] bufcopy;
|
|
|
|
// -----------------------------------------------------
|
|
// free local memory
|
|
// -----------------------------------------------------
|
|
|
|
memory->sfree(npartner);
|
|
memory->sfree(nshake);
|
|
memory->destroy_2d_int_array(partner_tag);
|
|
memory->destroy_2d_int_array(partner_mask);
|
|
memory->destroy_2d_int_array(partner_type);
|
|
memory->destroy_2d_int_array(partner_bondtype);
|
|
memory->destroy_2d_int_array(partner_shake);
|
|
memory->destroy_2d_int_array(partner_nshake);
|
|
|
|
// -----------------------------------------------------
|
|
// set bond_type and angle_type negative for SHAKE clusters
|
|
// must set for all SHAKE bonds and angles stored by each atom
|
|
// -----------------------------------------------------
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i] == 0) continue;
|
|
else if (shake_flag[i] == 1) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = anglefind(i,shake_atom[i][1],shake_atom[i][2]);
|
|
if (n >= 0) angle_type[i][n] = -angle_type[i][n];
|
|
} else if (shake_flag[i] == 2) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
} else if (shake_flag[i] == 3) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
} else if (shake_flag[i] == 4) {
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
n = bondfind(i,shake_atom[i][0],shake_atom[i][3]);
|
|
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
|
|
}
|
|
}
|
|
|
|
// -----------------------------------------------------
|
|
// print info on SHAKE clusters
|
|
// -----------------------------------------------------
|
|
|
|
int count1,count2,count3,count4;
|
|
count1 = count2 = count3 = count4 = 0;
|
|
for (i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i] == 1) count1++;
|
|
else if (shake_flag[i] == 2) count2++;
|
|
else if (shake_flag[i] == 3) count3++;
|
|
else if (shake_flag[i] == 4) count4++;
|
|
}
|
|
|
|
int tmp;
|
|
tmp = count1;
|
|
MPI_Allreduce(&tmp,&count1,1,MPI_INT,MPI_SUM,world);
|
|
tmp = count2;
|
|
MPI_Allreduce(&tmp,&count2,1,MPI_INT,MPI_SUM,world);
|
|
tmp = count3;
|
|
MPI_Allreduce(&tmp,&count3,1,MPI_INT,MPI_SUM,world);
|
|
tmp = count4;
|
|
MPI_Allreduce(&tmp,&count4,1,MPI_INT,MPI_SUM,world);
|
|
|
|
if (me == 0) {
|
|
if (screen) {
|
|
fprintf(screen," %d = # of size 2 clusters\n",count2/2);
|
|
fprintf(screen," %d = # of size 3 clusters\n",count3/3);
|
|
fprintf(screen," %d = # of size 4 clusters\n",count4/4);
|
|
fprintf(screen," %d = # of frozen angles\n",count1/3);
|
|
}
|
|
if (logfile) {
|
|
fprintf(logfile," %d = # of size 2 clusters\n",count2/2);
|
|
fprintf(logfile," %d = # of size 3 clusters\n",count3/3);
|
|
fprintf(logfile," %d = # of size 4 clusters\n",count4/4);
|
|
fprintf(logfile," %d = # of frozen angles\n",count1/3);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
update the unconstrained position of each atom
|
|
only for SHAKE clusters, else set to 0.0
|
|
assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::unconstrained_update()
|
|
{
|
|
double dtfmsq;
|
|
|
|
for (int i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i]) {
|
|
dtfmsq = dtfsq / mass[type[i]];
|
|
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
|
|
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
|
|
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
|
|
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
|
|
}
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void FixShake::shake2(int m)
|
|
{
|
|
// local atom IDs and constraint distances
|
|
|
|
int i0 = atom->map(shake_atom[m][0]);
|
|
int i1 = atom->map(shake_atom[m][1]);
|
|
double bond1 = bond_distance[shake_type[m][0]];
|
|
|
|
// r01 = distance vec between atoms, with PBC
|
|
|
|
double r01[3];
|
|
r01[0] = x[i0][0] - x[i1][0];
|
|
r01[1] = x[i0][1] - x[i1][1];
|
|
r01[2] = x[i0][2] - x[i1][2];
|
|
domain->minimum_image(r01);
|
|
|
|
// s01 = distance vec after unconstrained update, with PBC
|
|
|
|
double s01[3];
|
|
s01[0] = xshake[i0][0] - xshake[i1][0];
|
|
s01[1] = xshake[i0][1] - xshake[i1][1];
|
|
s01[2] = xshake[i0][2] - xshake[i1][2];
|
|
domain->minimum_image(s01);
|
|
|
|
// scalar distances between atoms
|
|
|
|
double r01sq = r01[0]*r01[0] + r01[1]*r01[1] + r01[2]*r01[2];
|
|
double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2];
|
|
|
|
// a,b,c = coeffs in quadratic equation for lamda
|
|
|
|
double invmass0 = 1.0/mass[type[i0]];
|
|
double invmass1 = 1.0/mass[type[i1]];
|
|
|
|
double a = (invmass0+invmass1)*(invmass0+invmass1) * r01sq;
|
|
double b = 2.0 * (invmass0+invmass1) *
|
|
(s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]);
|
|
double c = s01sq - bond1*bond1;
|
|
|
|
// error check
|
|
|
|
double determ = b*b - 4.0*a*c;
|
|
if (determ < 0.0) {
|
|
error->warning("Shake determinant < 0.0");
|
|
determ = 0.0;
|
|
}
|
|
|
|
// exact quadratic solution for lamda
|
|
|
|
double lamda,lamda1,lamda2;
|
|
lamda1 = (-b+sqrt(determ)) / (2.0*a);
|
|
lamda2 = (-b-sqrt(determ)) / (2.0*a);
|
|
|
|
if (fabs(lamda1) <= fabs(lamda2)) lamda = lamda1;
|
|
else lamda = lamda2;
|
|
|
|
// update forces if atom is owned by this processor
|
|
|
|
lamda /= dtfsq;
|
|
|
|
if (i0 < nlocal) {
|
|
f[i0][0] += lamda*r01[0];
|
|
f[i0][1] += lamda*r01[1];
|
|
f[i0][2] += lamda*r01[2];
|
|
}
|
|
|
|
if (i1 < nlocal) {
|
|
f[i1][0] -= lamda*r01[0];
|
|
f[i1][1] -= lamda*r01[1];
|
|
f[i1][2] -= lamda*r01[2];
|
|
}
|
|
|
|
if (vflag) {
|
|
int factor = 0;
|
|
if (i0 < nlocal) factor++;
|
|
if (i1 < nlocal) factor++;
|
|
double rfactor = 0.5 * factor;
|
|
|
|
virial[0] += rfactor*lamda*r01[0]*r01[0];
|
|
virial[1] += rfactor*lamda*r01[1]*r01[1];
|
|
virial[2] += rfactor*lamda*r01[2]*r01[2];
|
|
virial[3] += rfactor*lamda*r01[0]*r01[1];
|
|
virial[4] += rfactor*lamda*r01[0]*r01[2];
|
|
virial[5] += rfactor*lamda*r01[1]*r01[2];
|
|
}
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void FixShake::shake3(int m)
|
|
{
|
|
// 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]];
|
|
|
|
// r01,r02 = distance vec between atoms, with PBC
|
|
|
|
double r01[3];
|
|
r01[0] = x[i0][0] - x[i1][0];
|
|
r01[1] = x[i0][1] - x[i1][1];
|
|
r01[2] = x[i0][2] - x[i1][2];
|
|
domain->minimum_image(r01);
|
|
|
|
double r02[3];
|
|
r02[0] = x[i0][0] - x[i2][0];
|
|
r02[1] = x[i0][1] - x[i2][1];
|
|
r02[2] = x[i0][2] - x[i2][2];
|
|
domain->minimum_image(r02);
|
|
|
|
// s01,s02 = distance vec after unconstrained update, with PBC
|
|
|
|
double s01[3];
|
|
s01[0] = xshake[i0][0] - xshake[i1][0];
|
|
s01[1] = xshake[i0][1] - xshake[i1][1];
|
|
s01[2] = xshake[i0][2] - xshake[i1][2];
|
|
domain->minimum_image(s01);
|
|
|
|
double s02[3];
|
|
s02[0] = xshake[i0][0] - xshake[i2][0];
|
|
s02[1] = xshake[i0][1] - xshake[i2][1];
|
|
s02[2] = xshake[i0][2] - xshake[i2][2];
|
|
domain->minimum_image(s02);
|
|
|
|
// 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 s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2];
|
|
double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2];
|
|
|
|
// matrix coeffs and rhs for lamda equations
|
|
|
|
double invmass0 = 1.0/mass[type[i0]];
|
|
double invmass1 = 1.0/mass[type[i1]];
|
|
double invmass2 = 1.0/mass[type[i2]];
|
|
|
|
double a11 = 2.0 * (invmass0+invmass1) *
|
|
(s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]);
|
|
double a12 = 2.0 * invmass0 *
|
|
(s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]);
|
|
double a21 = 2.0 * invmass0 *
|
|
(s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]);
|
|
double a22 = 2.0 * (invmass0+invmass2) *
|
|
(s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]);
|
|
|
|
// inverse of matrix
|
|
|
|
double determ = a11*a22 - a12*a21;
|
|
if (determ == 0.0) error->one("Shake determinant = 0.0");
|
|
double determinv = 1.0/determ;
|
|
|
|
double a11inv = a22*determinv;
|
|
double a12inv = -a12*determinv;
|
|
double a21inv = -a21*determinv;
|
|
double a22inv = a11*determinv;
|
|
|
|
// quadratic correction coeffs
|
|
|
|
double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]);
|
|
|
|
double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq;
|
|
double quad1_0202 = invmass0*invmass0 * r02sq;
|
|
double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102;
|
|
|
|
double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq;
|
|
double quad2_0101 = invmass0*invmass0 * r01sq;
|
|
double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102;
|
|
|
|
// iterate until converged
|
|
|
|
double lamda01 = 0.0;
|
|
double lamda02 = 0.0;
|
|
int niter = 0;
|
|
int done = 0;
|
|
|
|
double quad1,quad2,b1,b2,lamda01_new,lamda02_new;
|
|
|
|
while (!done && niter < max_iter) {
|
|
quad1 = quad1_0101 * lamda01*lamda01 + quad1_0202 * lamda02*lamda02 +
|
|
quad1_0102 * lamda01*lamda02;
|
|
quad2 = quad2_0101 * lamda01*lamda01 + quad2_0202 * lamda02*lamda02 +
|
|
quad2_0102 * lamda01*lamda02;
|
|
|
|
b1 = bond1*bond1 - s01sq - quad1;
|
|
b2 = bond2*bond2 - s02sq - quad2;
|
|
|
|
lamda01_new = a11inv*b1 + a12inv*b2;
|
|
lamda02_new = a21inv*b1 + a22inv*b2;
|
|
|
|
done = 1;
|
|
if (fabs(lamda01_new-lamda01) > tolerance) done = 0;
|
|
if (fabs(lamda02_new-lamda02) > tolerance) done = 0;
|
|
|
|
lamda01 = lamda01_new;
|
|
lamda02 = lamda02_new;
|
|
niter++;
|
|
}
|
|
|
|
// update forces if atom is owned by this processor
|
|
|
|
lamda01 = lamda01/dtfsq;
|
|
lamda02 = lamda02/dtfsq;
|
|
|
|
if (i0 < nlocal) {
|
|
f[i0][0] += lamda01*r01[0] + lamda02*r02[0];
|
|
f[i0][1] += lamda01*r01[1] + lamda02*r02[1];
|
|
f[i0][2] += lamda01*r01[2] + lamda02*r02[2];
|
|
}
|
|
|
|
if (i1 < nlocal) {
|
|
f[i1][0] -= lamda01*r01[0];
|
|
f[i1][1] -= lamda01*r01[1];
|
|
f[i1][2] -= lamda01*r01[2];
|
|
}
|
|
|
|
if (i2 < nlocal) {
|
|
f[i2][0] -= lamda02*r02[0];
|
|
f[i2][1] -= lamda02*r02[1];
|
|
f[i2][2] -= lamda02*r02[2];
|
|
}
|
|
|
|
if (vflag) {
|
|
int factor = 0;
|
|
if (i0 < nlocal) factor++;
|
|
if (i1 < nlocal) factor++;
|
|
if (i2 < nlocal) factor++;
|
|
double rfactor = factor/3.0;
|
|
|
|
virial[0] += rfactor*lamda01*r01[0]*r01[0];
|
|
virial[1] += rfactor*lamda01*r01[1]*r01[1];
|
|
virial[2] += rfactor*lamda01*r01[2]*r01[2];
|
|
virial[3] += rfactor*lamda01*r01[0]*r01[1];
|
|
virial[4] += rfactor*lamda01*r01[0]*r01[2];
|
|
virial[5] += rfactor*lamda01*r01[1]*r01[2];
|
|
|
|
virial[0] += rfactor*lamda02*r02[0]*r02[0];
|
|
virial[1] += rfactor*lamda02*r02[1]*r02[1];
|
|
virial[2] += rfactor*lamda02*r02[2]*r02[2];
|
|
virial[3] += rfactor*lamda02*r02[0]*r02[1];
|
|
virial[4] += rfactor*lamda02*r02[0]*r02[2];
|
|
virial[5] += rfactor*lamda02*r02[1]*r02[2];
|
|
}
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void FixShake::shake4(int m)
|
|
{
|
|
// 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]);
|
|
int i3 = atom->map(shake_atom[m][3]);
|
|
double bond1 = bond_distance[shake_type[m][0]];
|
|
double bond2 = bond_distance[shake_type[m][1]];
|
|
double bond3 = bond_distance[shake_type[m][2]];
|
|
|
|
// r01,r02,r03 = distance vec between atoms, with PBC
|
|
|
|
double r01[3];
|
|
r01[0] = x[i0][0] - x[i1][0];
|
|
r01[1] = x[i0][1] - x[i1][1];
|
|
r01[2] = x[i0][2] - x[i1][2];
|
|
domain->minimum_image(r01);
|
|
|
|
double r02[3];
|
|
r02[0] = x[i0][0] - x[i2][0];
|
|
r02[1] = x[i0][1] - x[i2][1];
|
|
r02[2] = x[i0][2] - x[i2][2];
|
|
domain->minimum_image(r02);
|
|
|
|
double r03[3];
|
|
r03[0] = x[i0][0] - x[i3][0];
|
|
r03[1] = x[i0][1] - x[i3][1];
|
|
r03[2] = x[i0][2] - x[i3][2];
|
|
domain->minimum_image(r03);
|
|
|
|
// s01,s02,s03 = distance vec after unconstrained update, with PBC
|
|
|
|
double s01[3];
|
|
s01[0] = xshake[i0][0] - xshake[i1][0];
|
|
s01[1] = xshake[i0][1] - xshake[i1][1];
|
|
s01[2] = xshake[i0][2] - xshake[i1][2];
|
|
domain->minimum_image(s01);
|
|
|
|
double s02[3];
|
|
s02[0] = xshake[i0][0] - xshake[i2][0];
|
|
s02[1] = xshake[i0][1] - xshake[i2][1];
|
|
s02[2] = xshake[i0][2] - xshake[i2][2];
|
|
domain->minimum_image(s02);
|
|
|
|
double s03[3];
|
|
s03[0] = xshake[i0][0] - xshake[i3][0];
|
|
s03[1] = xshake[i0][1] - xshake[i3][1];
|
|
s03[2] = xshake[i0][2] - xshake[i3][2];
|
|
domain->minimum_image(s03);
|
|
|
|
// 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 r03sq = r03[0]*r03[0] + r03[1]*r03[1] + r03[2]*r03[2];
|
|
double s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2];
|
|
double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2];
|
|
double s03sq = s03[0]*s03[0] + s03[1]*s03[1] + s03[2]*s03[2];
|
|
|
|
// matrix coeffs and rhs for lamda equations
|
|
|
|
double invmass0 = 1.0/mass[type[i0]];
|
|
double invmass1 = 1.0/mass[type[i1]];
|
|
double invmass2 = 1.0/mass[type[i2]];
|
|
double invmass3 = 1.0/mass[type[i3]];
|
|
|
|
double a11 = 2.0 * (invmass0+invmass1) *
|
|
(s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]);
|
|
double a12 = 2.0 * invmass0 *
|
|
(s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]);
|
|
double a13 = 2.0 * invmass0 *
|
|
(s01[0]*r03[0] + s01[1]*r03[1] + s01[2]*r03[2]);
|
|
double a21 = 2.0 * invmass0 *
|
|
(s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]);
|
|
double a22 = 2.0 * (invmass0+invmass2) *
|
|
(s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]);
|
|
double a23 = 2.0 * invmass0 *
|
|
(s02[0]*r03[0] + s02[1]*r03[1] + s02[2]*r03[2]);
|
|
double a31 = 2.0 * invmass0 *
|
|
(s03[0]*r01[0] + s03[1]*r01[1] + s03[2]*r01[2]);
|
|
double a32 = 2.0 * invmass0 *
|
|
(s03[0]*r02[0] + s03[1]*r02[1] + s03[2]*r02[2]);
|
|
double a33 = 2.0 * (invmass0+invmass3) *
|
|
(s03[0]*r03[0] + s03[1]*r03[1] + s03[2]*r03[2]);
|
|
|
|
// inverse of matrix;
|
|
|
|
double determ = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 -
|
|
a11*a23*a32 - a12*a21*a33 - a13*a22*a31;
|
|
if (determ == 0.0) error->one("Shake determinant = 0.0");
|
|
double determinv = 1.0/determ;
|
|
|
|
double a11inv = determinv * (a22*a33 - a23*a32);
|
|
double a12inv = -determinv * (a12*a33 - a13*a32);
|
|
double a13inv = determinv * (a12*a23 - a13*a22);
|
|
double a21inv = -determinv * (a21*a33 - a23*a31);
|
|
double a22inv = determinv * (a11*a33 - a13*a31);
|
|
double a23inv = -determinv * (a11*a23 - a13*a21);
|
|
double a31inv = determinv * (a21*a32 - a22*a31);
|
|
double a32inv = -determinv * (a11*a32 - a12*a31);
|
|
double a33inv = determinv * (a11*a22 - a12*a21);
|
|
|
|
// quadratic correction coeffs
|
|
|
|
double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]);
|
|
double r0103 = (r01[0]*r03[0] + r01[1]*r03[1] + r01[2]*r03[2]);
|
|
double r0203 = (r02[0]*r03[0] + r02[1]*r03[1] + r02[2]*r03[2]);
|
|
|
|
double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq;
|
|
double quad1_0202 = invmass0*invmass0 * r02sq;
|
|
double quad1_0303 = invmass0*invmass0 * r03sq;
|
|
double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102;
|
|
double quad1_0103 = 2.0 * (invmass0+invmass1)*invmass0 * r0103;
|
|
double quad1_0203 = 2.0 * invmass0*invmass0 * r0203;
|
|
|
|
double quad2_0101 = invmass0*invmass0 * r01sq;
|
|
double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq;
|
|
double quad2_0303 = invmass0*invmass0 * r03sq;
|
|
double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102;
|
|
double quad2_0103 = 2.0 * invmass0*invmass0 * r0103;
|
|
double quad2_0203 = 2.0 * (invmass0+invmass2)*invmass0 * r0203;
|
|
|
|
double quad3_0101 = invmass0*invmass0 * r01sq;
|
|
double quad3_0202 = invmass0*invmass0 * r02sq;
|
|
double quad3_0303 = (invmass0+invmass3)*(invmass0+invmass3) * r03sq;
|
|
double quad3_0102 = 2.0 * invmass0*invmass0 * r0102;
|
|
double quad3_0103 = 2.0 * (invmass0+invmass3)*invmass0 * r0103;
|
|
double quad3_0203 = 2.0 * (invmass0+invmass3)*invmass0 * r0203;
|
|
|
|
// iterate until converged
|
|
|
|
double lamda01 = 0.0;
|
|
double lamda02 = 0.0;
|
|
double lamda03 = 0.0;
|
|
int niter = 0;
|
|
int done = 0;
|
|
|
|
double quad1,quad2,quad3,b1,b2,b3,lamda01_new,lamda02_new,lamda03_new;
|
|
|
|
while (!done && niter < max_iter) {
|
|
quad1 = quad1_0101 * lamda01*lamda01 +
|
|
quad1_0202 * lamda02*lamda02 +
|
|
quad1_0303 * lamda03*lamda03 +
|
|
quad1_0102 * lamda01*lamda02 +
|
|
quad1_0103 * lamda01*lamda03 +
|
|
quad1_0203 * lamda02*lamda03;
|
|
|
|
quad2 = quad2_0101 * lamda01*lamda01 +
|
|
quad2_0202 * lamda02*lamda02 +
|
|
quad2_0303 * lamda03*lamda03 +
|
|
quad2_0102 * lamda01*lamda02 +
|
|
quad2_0103 * lamda01*lamda03 +
|
|
quad2_0203 * lamda02*lamda03;
|
|
|
|
quad3 = quad3_0101 * lamda01*lamda01 +
|
|
quad3_0202 * lamda02*lamda02 +
|
|
quad3_0303 * lamda03*lamda03 +
|
|
quad3_0102 * lamda01*lamda02 +
|
|
quad3_0103 * lamda01*lamda03 +
|
|
quad3_0203 * lamda02*lamda03;
|
|
|
|
b1 = bond1*bond1 - s01sq - quad1;
|
|
b2 = bond2*bond2 - s02sq - quad2;
|
|
b3 = bond3*bond3 - s03sq - quad3;
|
|
|
|
lamda01_new = a11inv*b1 + a12inv*b2 + a13inv*b3;
|
|
lamda02_new = a21inv*b1 + a22inv*b2 + a23inv*b3;
|
|
lamda03_new = a31inv*b1 + a32inv*b2 + a33inv*b3;
|
|
|
|
done = 1;
|
|
if (fabs(lamda01_new-lamda01) > tolerance) done = 0;
|
|
if (fabs(lamda02_new-lamda02) > tolerance) done = 0;
|
|
if (fabs(lamda03_new-lamda03) > tolerance) done = 0;
|
|
|
|
lamda01 = lamda01_new;
|
|
lamda02 = lamda02_new;
|
|
lamda03 = lamda03_new;
|
|
niter++;
|
|
}
|
|
|
|
// update forces if atom is owned by this processor
|
|
|
|
lamda01 = lamda01/dtfsq;
|
|
lamda02 = lamda02/dtfsq;
|
|
lamda03 = lamda03/dtfsq;
|
|
|
|
if (i0 < nlocal) {
|
|
f[i0][0] += lamda01*r01[0] + lamda02*r02[0] + lamda03*r03[0];
|
|
f[i0][1] += lamda01*r01[1] + lamda02*r02[1] + lamda03*r03[1];
|
|
f[i0][2] += lamda01*r01[2] + lamda02*r02[2] + lamda03*r03[2];
|
|
}
|
|
|
|
if (i1 < nlocal) {
|
|
f[i1][0] -= lamda01*r01[0];
|
|
f[i1][1] -= lamda01*r01[1];
|
|
f[i1][2] -= lamda01*r01[2];
|
|
}
|
|
|
|
if (i2 < nlocal) {
|
|
f[i2][0] -= lamda02*r02[0];
|
|
f[i2][1] -= lamda02*r02[1];
|
|
f[i2][2] -= lamda02*r02[2];
|
|
}
|
|
|
|
if (i3 < nlocal) {
|
|
f[i3][0] -= lamda03*r03[0];
|
|
f[i3][1] -= lamda03*r03[1];
|
|
f[i3][2] -= lamda03*r03[2];
|
|
}
|
|
|
|
if (vflag) {
|
|
int factor = 0;
|
|
if (i0 < nlocal) factor++;
|
|
if (i1 < nlocal) factor++;
|
|
if (i2 < nlocal) factor++;
|
|
if (i3 < nlocal) factor++;
|
|
double rfactor = 0.25*factor;
|
|
|
|
virial[0] += rfactor*lamda01*r01[0]*r01[0];
|
|
virial[1] += rfactor*lamda01*r01[1]*r01[1];
|
|
virial[2] += rfactor*lamda01*r01[2]*r01[2];
|
|
virial[3] += rfactor*lamda01*r01[0]*r01[1];
|
|
virial[4] += rfactor*lamda01*r01[0]*r01[2];
|
|
virial[5] += rfactor*lamda01*r01[1]*r01[2];
|
|
|
|
virial[0] += rfactor*lamda02*r02[0]*r02[0];
|
|
virial[1] += rfactor*lamda02*r02[1]*r02[1];
|
|
virial[2] += rfactor*lamda02*r02[2]*r02[2];
|
|
virial[3] += rfactor*lamda02*r02[0]*r02[1];
|
|
virial[4] += rfactor*lamda02*r02[0]*r02[2];
|
|
virial[5] += rfactor*lamda02*r02[1]*r02[2];
|
|
|
|
virial[0] += rfactor*lamda03*r03[0]*r03[0];
|
|
virial[1] += rfactor*lamda03*r03[1]*r03[1];
|
|
virial[2] += rfactor*lamda03*r03[2]*r03[2];
|
|
virial[3] += rfactor*lamda03*r03[0]*r03[1];
|
|
virial[4] += rfactor*lamda03*r03[0]*r03[2];
|
|
virial[5] += rfactor*lamda03*r03[1]*r03[2];
|
|
}
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void FixShake::shake3angle(int m)
|
|
{
|
|
// 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[i0][0] - x[i1][0];
|
|
r01[1] = x[i0][1] - x[i1][1];
|
|
r01[2] = x[i0][2] - x[i1][2];
|
|
domain->minimum_image(r01);
|
|
|
|
double r02[3];
|
|
r02[0] = x[i0][0] - x[i2][0];
|
|
r02[1] = x[i0][1] - x[i2][1];
|
|
r02[2] = x[i0][2] - x[i2][2];
|
|
domain->minimum_image(r02);
|
|
|
|
double r12[3];
|
|
r12[0] = x[i1][0] - x[i2][0];
|
|
r12[1] = x[i1][1] - x[i2][1];
|
|
r12[2] = x[i1][2] - x[i2][2];
|
|
domain->minimum_image(r12);
|
|
|
|
// s01,s02,s12 = distance vec after unconstrained update, with PBC
|
|
|
|
double s01[3];
|
|
s01[0] = xshake[i0][0] - xshake[i1][0];
|
|
s01[1] = xshake[i0][1] - xshake[i1][1];
|
|
s01[2] = xshake[i0][2] - xshake[i1][2];
|
|
domain->minimum_image(s01);
|
|
|
|
double s02[3];
|
|
s02[0] = xshake[i0][0] - xshake[i2][0];
|
|
s02[1] = xshake[i0][1] - xshake[i2][1];
|
|
s02[2] = xshake[i0][2] - xshake[i2][2];
|
|
domain->minimum_image(s02);
|
|
|
|
double s12[3];
|
|
s12[0] = xshake[i1][0] - xshake[i2][0];
|
|
s12[1] = xshake[i1][1] - xshake[i2][1];
|
|
s12[2] = xshake[i1][2] - xshake[i2][2];
|
|
domain->minimum_image(s12);
|
|
|
|
// 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 s01sq = s01[0]*s01[0] + s01[1]*s01[1] + s01[2]*s01[2];
|
|
double s02sq = s02[0]*s02[0] + s02[1]*s02[1] + s02[2]*s02[2];
|
|
double s12sq = s12[0]*s12[0] + s12[1]*s12[1] + s12[2]*s12[2];
|
|
|
|
// matrix coeffs and rhs for lamda equations
|
|
|
|
double invmass0 = 1.0/mass[type[i0]];
|
|
double invmass1 = 1.0/mass[type[i1]];
|
|
double invmass2 = 1.0/mass[type[i2]];
|
|
|
|
double a11 = 2.0 * (invmass0+invmass1) *
|
|
(s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]);
|
|
double a12 = 2.0 * invmass0 *
|
|
(s01[0]*r02[0] + s01[1]*r02[1] + s01[2]*r02[2]);
|
|
double a13 = - 2.0 * invmass1 *
|
|
(s01[0]*r12[0] + s01[1]*r12[1] + s01[2]*r12[2]);
|
|
double a21 = 2.0 * invmass0 *
|
|
(s02[0]*r01[0] + s02[1]*r01[1] + s02[2]*r01[2]);
|
|
double a22 = 2.0 * (invmass0+invmass2) *
|
|
(s02[0]*r02[0] + s02[1]*r02[1] + s02[2]*r02[2]);
|
|
double a23 = 2.0 * invmass2 *
|
|
(s02[0]*r12[0] + s02[1]*r12[1] + s02[2]*r12[2]);
|
|
double a31 = - 2.0 * invmass1 *
|
|
(s12[0]*r01[0] + s12[1]*r01[1] + s12[2]*r01[2]);
|
|
double a32 = 2.0 * invmass2 *
|
|
(s12[0]*r02[0] + s12[1]*r02[1] + s12[2]*r02[2]);
|
|
double a33 = 2.0 * (invmass1+invmass2) *
|
|
(s12[0]*r12[0] + s12[1]*r12[1] + s12[2]*r12[2]);
|
|
|
|
// inverse of matrix
|
|
|
|
double determ = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 -
|
|
a11*a23*a32 - a12*a21*a33 - a13*a22*a31;
|
|
if (determ == 0.0) error->one("Shake determinant = 0.0");
|
|
double determinv = 1.0/determ;
|
|
|
|
double a11inv = determinv * (a22*a33 - a23*a32);
|
|
double a12inv = -determinv * (a12*a33 - a13*a32);
|
|
double a13inv = determinv * (a12*a23 - a13*a22);
|
|
double a21inv = -determinv * (a21*a33 - a23*a31);
|
|
double a22inv = determinv * (a11*a33 - a13*a31);
|
|
double a23inv = -determinv * (a11*a23 - a13*a21);
|
|
double a31inv = determinv * (a21*a32 - a22*a31);
|
|
double a32inv = -determinv * (a11*a32 - a12*a31);
|
|
double a33inv = determinv * (a11*a22 - a12*a21);
|
|
|
|
// quadratic correction coeffs
|
|
|
|
double r0102 = (r01[0]*r02[0] + r01[1]*r02[1] + r01[2]*r02[2]);
|
|
double r0112 = (r01[0]*r12[0] + r01[1]*r12[1] + r01[2]*r12[2]);
|
|
double r0212 = (r02[0]*r12[0] + r02[1]*r12[1] + r02[2]*r12[2]);
|
|
|
|
double quad1_0101 = (invmass0+invmass1)*(invmass0+invmass1) * r01sq;
|
|
double quad1_0202 = invmass0*invmass0 * r02sq;
|
|
double quad1_1212 = invmass1*invmass1 * r12sq;
|
|
double quad1_0102 = 2.0 * (invmass0+invmass1)*invmass0 * r0102;
|
|
double quad1_0112 = - 2.0 * (invmass0+invmass1)*invmass1 * r0112;
|
|
double quad1_0212 = - 2.0 * invmass0*invmass1 * r0212;
|
|
|
|
double quad2_0101 = invmass0*invmass0 * r01sq;
|
|
double quad2_0202 = (invmass0+invmass2)*(invmass0+invmass2) * r02sq;
|
|
double quad2_1212 = invmass2*invmass2 * r12sq;
|
|
double quad2_0102 = 2.0 * (invmass0+invmass2)*invmass0 * r0102;
|
|
double quad2_0112 = 2.0 * invmass0*invmass2 * r0112;
|
|
double quad2_0212 = 2.0 * (invmass0+invmass2)*invmass2 * r0212;
|
|
|
|
double quad3_0101 = invmass1*invmass1 * r01sq;
|
|
double quad3_0202 = invmass2*invmass2 * r02sq;
|
|
double quad3_1212 = (invmass1+invmass2)*(invmass1+invmass2) * r12sq;
|
|
double quad3_0102 = - 2.0 * invmass1*invmass2 * r0102;
|
|
double quad3_0112 = - 2.0 * (invmass1+invmass2)*invmass1 * r0112;
|
|
double quad3_0212 = 2.0 * (invmass1+invmass2)*invmass2 * r0212;
|
|
|
|
// iterate until converged
|
|
|
|
double lamda01 = 0.0;
|
|
double lamda02 = 0.0;
|
|
double lamda12 = 0.0;
|
|
int niter = 0;
|
|
int done = 0;
|
|
|
|
double quad1,quad2,quad3,b1,b2,b3,lamda01_new,lamda02_new,lamda12_new;
|
|
|
|
while (!done && niter < max_iter) {
|
|
quad1 = quad1_0101 * lamda01*lamda01 +
|
|
quad1_0202 * lamda02*lamda02 +
|
|
quad1_1212 * lamda12*lamda12 +
|
|
quad1_0102 * lamda01*lamda02 +
|
|
quad1_0112 * lamda01*lamda12 +
|
|
quad1_0212 * lamda02*lamda12;
|
|
|
|
quad2 = quad2_0101 * lamda01*lamda01 +
|
|
quad2_0202 * lamda02*lamda02 +
|
|
quad2_1212 * lamda12*lamda12 +
|
|
quad2_0102 * lamda01*lamda02 +
|
|
quad2_0112 * lamda01*lamda12 +
|
|
quad2_0212 * lamda02*lamda12;
|
|
|
|
quad3 = quad3_0101 * lamda01*lamda01 +
|
|
quad3_0202 * lamda02*lamda02 +
|
|
quad3_1212 * lamda12*lamda12 +
|
|
quad3_0102 * lamda01*lamda02 +
|
|
quad3_0112 * lamda01*lamda12 +
|
|
quad3_0212 * lamda02*lamda12;
|
|
|
|
b1 = bond1*bond1 - s01sq - quad1;
|
|
b2 = bond2*bond2 - s02sq - quad2;
|
|
b3 = bond12*bond12 - s12sq - quad3;
|
|
|
|
lamda01_new = a11inv*b1 + a12inv*b2 + a13inv*b3;
|
|
lamda02_new = a21inv*b1 + a22inv*b2 + a23inv*b3;
|
|
lamda12_new = a31inv*b1 + a32inv*b2 + a33inv*b3;
|
|
|
|
done = 1;
|
|
if (fabs(lamda01_new-lamda01) > tolerance) done = 0;
|
|
if (fabs(lamda02_new-lamda02) > tolerance) done = 0;
|
|
if (fabs(lamda12_new-lamda12) > tolerance) done = 0;
|
|
|
|
lamda01 = lamda01_new;
|
|
lamda02 = lamda02_new;
|
|
lamda12 = lamda12_new;
|
|
niter++;
|
|
}
|
|
|
|
// update forces if atom is owned by this processor
|
|
|
|
lamda01 = lamda01/dtfsq;
|
|
lamda02 = lamda02/dtfsq;
|
|
lamda12 = lamda12/dtfsq;
|
|
|
|
if (i0 < nlocal) {
|
|
f[i0][0] += lamda01*r01[0] + lamda02*r02[0];
|
|
f[i0][1] += lamda01*r01[1] + lamda02*r02[1];
|
|
f[i0][2] += lamda01*r01[2] + lamda02*r02[2];
|
|
}
|
|
|
|
if (i1 < nlocal) {
|
|
f[i1][0] -= lamda01*r01[0] - lamda12*r12[0];
|
|
f[i1][1] -= lamda01*r01[1] - lamda12*r12[1];
|
|
f[i1][2] -= lamda01*r01[2] - lamda12*r12[2];
|
|
}
|
|
|
|
if (i2 < nlocal) {
|
|
f[i2][0] -= lamda02*r02[0] + lamda12*r12[0];
|
|
f[i2][1] -= lamda02*r02[1] + lamda12*r12[1];
|
|
f[i2][2] -= lamda02*r02[2] + lamda12*r12[2];
|
|
}
|
|
|
|
if (vflag) {
|
|
int factor = 0;
|
|
if (i0 < nlocal) factor++;
|
|
if (i1 < nlocal) factor++;
|
|
if (i2 < nlocal) factor++;
|
|
double rfactor = factor/3.0;
|
|
|
|
virial[0] += rfactor*lamda01*r01[0]*r01[0];
|
|
virial[1] += rfactor*lamda01*r01[1]*r01[1];
|
|
virial[2] += rfactor*lamda01*r01[2]*r01[2];
|
|
virial[3] += rfactor*lamda01*r01[0]*r01[1];
|
|
virial[4] += rfactor*lamda01*r01[0]*r01[2];
|
|
virial[5] += rfactor*lamda01*r01[1]*r01[2];
|
|
|
|
virial[0] += rfactor*lamda02*r02[0]*r02[0];
|
|
virial[1] += rfactor*lamda02*r02[1]*r02[1];
|
|
virial[2] += rfactor*lamda02*r02[2]*r02[2];
|
|
virial[3] += rfactor*lamda02*r02[0]*r02[1];
|
|
virial[4] += rfactor*lamda02*r02[0]*r02[2];
|
|
virial[5] += rfactor*lamda02*r02[1]*r02[2];
|
|
|
|
virial[0] += rfactor*lamda12*r12[0]*r12[0];
|
|
virial[1] += rfactor*lamda12*r12[1]*r12[1];
|
|
virial[2] += rfactor*lamda12*r12[2]*r12[2];
|
|
virial[3] += rfactor*lamda12*r12[0]*r12[1];
|
|
virial[4] += rfactor*lamda12*r12[0]*r12[2];
|
|
virial[5] += rfactor*lamda12*r12[1]*r12[2];
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
print-out bond & angle statistics
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::stats()
|
|
{
|
|
int i,j,m,n,iatom,jatom,katom;
|
|
double delx,dely,delz;
|
|
double r,r1,r2,r3,angle;
|
|
|
|
// zero out accumulators
|
|
|
|
int nb = atom->nbondtypes + 1;
|
|
int na = atom->nangletypes + 1;
|
|
|
|
for (i = 0; i < nb; i++) {
|
|
b_count[i] = 0;
|
|
b_ave[i] = b_max[i] = 0.0;
|
|
b_min[i] = BIG;
|
|
}
|
|
for (i = 0; i < na; i++) {
|
|
a_count[i] = 0;
|
|
a_ave[i] = a_max[i] = 0.0;
|
|
a_min[i] = BIG;
|
|
}
|
|
|
|
// log stats for each bond & angle
|
|
// OK to double count since are just averaging
|
|
|
|
double **x = atom->x;
|
|
int nlocal = atom->nlocal;
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i] == 0) continue;
|
|
|
|
// bond stats
|
|
|
|
n = shake_flag[i];
|
|
if (n == 1) n = 3;
|
|
iatom = atom->map(shake_atom[i][0]);
|
|
for (j = 1; j < n; j++) {
|
|
jatom = atom->map(shake_atom[i][j]);
|
|
delx = x[iatom][0] - x[jatom][0];
|
|
dely = x[iatom][1] - x[jatom][1];
|
|
delz = x[iatom][2] - x[jatom][2];
|
|
domain->minimum_image(&delx,&dely,&delz);
|
|
r = sqrt(delx*delx + dely*dely + delz*delz);
|
|
|
|
m = shake_type[i][j-1];
|
|
b_count[m]++;
|
|
b_ave[m] += r;
|
|
b_max[m] = MAX(b_max[m],r);
|
|
b_min[m] = MIN(b_min[m],r);
|
|
}
|
|
|
|
// angle stats
|
|
|
|
if (shake_flag[i] == 1) {
|
|
iatom = atom->map(shake_atom[i][0]);
|
|
jatom = atom->map(shake_atom[i][1]);
|
|
katom = atom->map(shake_atom[i][2]);
|
|
|
|
delx = x[iatom][0] - x[jatom][0];
|
|
dely = x[iatom][1] - x[jatom][1];
|
|
delz = x[iatom][2] - x[jatom][2];
|
|
domain->minimum_image(&delx,&dely,&delz);
|
|
r1 = sqrt(delx*delx + dely*dely + delz*delz);
|
|
|
|
delx = x[iatom][0] - x[katom][0];
|
|
dely = x[iatom][1] - x[katom][1];
|
|
delz = x[iatom][2] - x[katom][2];
|
|
domain->minimum_image(&delx,&dely,&delz);
|
|
r2 = sqrt(delx*delx + dely*dely + delz*delz);
|
|
|
|
delx = x[jatom][0] - x[katom][0];
|
|
dely = x[jatom][1] - x[katom][1];
|
|
delz = x[jatom][2] - x[katom][2];
|
|
domain->minimum_image(&delx,&dely,&delz);
|
|
r3 = sqrt(delx*delx + dely*dely + delz*delz);
|
|
|
|
angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2));
|
|
angle *= 180.0/PI;
|
|
m = shake_type[i][2];
|
|
a_count[m]++;
|
|
a_ave[m] += angle;
|
|
a_max[m] = MAX(a_max[m],angle);
|
|
a_min[m] = MIN(a_min[m],angle);
|
|
}
|
|
}
|
|
|
|
// sum across all procs
|
|
|
|
MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world);
|
|
MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world);
|
|
MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world);
|
|
MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world);
|
|
|
|
MPI_Allreduce(a_count,a_count_all,na,MPI_INT,MPI_SUM,world);
|
|
MPI_Allreduce(a_ave,a_ave_all,na,MPI_DOUBLE,MPI_SUM,world);
|
|
MPI_Allreduce(a_max,a_max_all,na,MPI_DOUBLE,MPI_MAX,world);
|
|
MPI_Allreduce(a_min,a_min_all,na,MPI_DOUBLE,MPI_MIN,world);
|
|
|
|
// print stats only for non-zero counts
|
|
|
|
if (me == 0) {
|
|
if (screen) {
|
|
fprintf(screen,"SHAKE stats (type/ave/delta) on step %d\n",
|
|
update->ntimestep);
|
|
for (i = 1; i < nb; i++)
|
|
if (b_count_all[i])
|
|
fprintf(screen," %d %g %g\n",i,
|
|
b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i]);
|
|
for (i = 1; i < na; i++)
|
|
if (a_count_all[i])
|
|
fprintf(screen," %d %g %g\n",i,
|
|
a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]);
|
|
}
|
|
if (logfile) {
|
|
fprintf(logfile,"SHAKE stats (type/ave/delta) on step %d\n",
|
|
update->ntimestep);
|
|
for (i = 0; i < nb; i++)
|
|
if (b_count_all[i])
|
|
fprintf(logfile," %d %g %g\n",i,
|
|
b_ave_all[i]/b_count_all[i],b_max_all[i]-b_min_all[i]);
|
|
for (i = 0; i < na; i++)
|
|
if (a_count_all[i])
|
|
fprintf(logfile," %d %g %g\n",i,
|
|
a_ave_all[i]/a_count_all[i],a_max_all[i]-a_min_all[i]);
|
|
}
|
|
}
|
|
|
|
// next timestep for stats
|
|
|
|
next_output += output_every;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
find a bond between global tags n1 and n2 stored with local atom i
|
|
return -1 if don't find it
|
|
return bond index if do find it
|
|
------------------------------------------------------------------------- */
|
|
|
|
int FixShake::bondfind(int i, int n1, int n2)
|
|
{
|
|
int *tag = atom->tag;
|
|
int **bond_atom = atom->bond_atom;
|
|
int nbonds = atom->num_bond[i];
|
|
|
|
int m;
|
|
for (m = 0; m < nbonds; m++) {
|
|
if (n1 == tag[i] && n2 == bond_atom[i][m]) break;
|
|
if (n1 == bond_atom[i][m] && n2 == tag[i]) break;
|
|
}
|
|
if (m < nbonds) return m;
|
|
return -1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
find an angle with global end atoms n1 and n2 stored with local atom i
|
|
return -1 if don't find it
|
|
return angle index if do find it
|
|
------------------------------------------------------------------------- */
|
|
|
|
int FixShake::anglefind(int i, int n1, int n2)
|
|
{
|
|
int **angle_atom1 = atom->angle_atom1;
|
|
int **angle_atom3 = atom->angle_atom3;
|
|
int nangles = atom->num_angle[i];
|
|
|
|
int m;
|
|
for (m = 0; m < nangles; m++) {
|
|
if (n1 == angle_atom1[i][m] && n2 == angle_atom3[i][m]) break;
|
|
if (n1 == angle_atom3[i][m] && n2 == angle_atom1[i][m]) break;
|
|
}
|
|
if (m < nangles) return m;
|
|
return -1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
memory usage of local atom-based arrays
|
|
------------------------------------------------------------------------- */
|
|
|
|
int FixShake::memory_usage()
|
|
{
|
|
int nmax = atom->nmax;
|
|
int bytes = nmax * sizeof(int);
|
|
bytes += nmax*4 * sizeof(int);
|
|
bytes += nmax*3 * sizeof(int);
|
|
bytes += nmax*3 * sizeof(double);
|
|
return bytes;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
allocate local atom-based arrays
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::grow_arrays(int nmax)
|
|
{
|
|
shake_flag = (int *)
|
|
memory->srealloc(shake_flag,nmax*sizeof(int),"shake:shake_flag");
|
|
shake_atom =
|
|
memory->grow_2d_int_array(shake_atom,nmax,4,"shake:shake_atom");
|
|
shake_type =
|
|
memory->grow_2d_int_array(shake_type,nmax,3,"shake:shake_type");
|
|
memory->destroy_2d_double_array(xshake);
|
|
xshake = memory->create_2d_double_array(nmax,3,"shake:xshake");
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
copy values within local atom-based arrays
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::copy_arrays(int i, int j)
|
|
{
|
|
int flag = shake_flag[j] = shake_flag[i];
|
|
if (flag == 1) {
|
|
shake_atom[j][0] = shake_atom[i][0];
|
|
shake_atom[j][1] = shake_atom[i][1];
|
|
shake_atom[j][2] = shake_atom[i][2];
|
|
shake_type[j][0] = shake_type[i][0];
|
|
shake_type[j][1] = shake_type[i][1];
|
|
shake_type[j][2] = shake_type[i][2];
|
|
} else if (flag == 2) {
|
|
shake_atom[j][0] = shake_atom[i][0];
|
|
shake_atom[j][1] = shake_atom[i][1];
|
|
shake_type[j][0] = shake_type[i][0];
|
|
} else if (flag == 3) {
|
|
shake_atom[j][0] = shake_atom[i][0];
|
|
shake_atom[j][1] = shake_atom[i][1];
|
|
shake_atom[j][2] = shake_atom[i][2];
|
|
shake_type[j][0] = shake_type[i][0];
|
|
shake_type[j][1] = shake_type[i][1];
|
|
} else if (flag == 4) {
|
|
shake_atom[j][0] = shake_atom[i][0];
|
|
shake_atom[j][1] = shake_atom[i][1];
|
|
shake_atom[j][2] = shake_atom[i][2];
|
|
shake_atom[j][3] = shake_atom[i][3];
|
|
shake_type[j][0] = shake_type[i][0];
|
|
shake_type[j][1] = shake_type[i][1];
|
|
shake_type[j][2] = shake_type[i][2];
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
pack values in local atom-based arrays for exchange with another proc
|
|
------------------------------------------------------------------------- */
|
|
|
|
int FixShake::pack_exchange(int i, double *buf)
|
|
{
|
|
int m = 0;
|
|
buf[m++] = shake_flag[i];
|
|
int flag = shake_flag[i];
|
|
if (flag == 1) {
|
|
buf[m++] = shake_atom[i][0];
|
|
buf[m++] = shake_atom[i][1];
|
|
buf[m++] = shake_atom[i][2];
|
|
buf[m++] = shake_type[i][0];
|
|
buf[m++] = shake_type[i][1];
|
|
buf[m++] = shake_type[i][2];
|
|
} else if (flag == 2) {
|
|
buf[m++] = shake_atom[i][0];
|
|
buf[m++] = shake_atom[i][1];
|
|
buf[m++] = shake_type[i][0];
|
|
} else if (flag == 3) {
|
|
buf[m++] = shake_atom[i][0];
|
|
buf[m++] = shake_atom[i][1];
|
|
buf[m++] = shake_atom[i][2];
|
|
buf[m++] = shake_type[i][0];
|
|
buf[m++] = shake_type[i][1];
|
|
} else if (flag == 4) {
|
|
buf[m++] = shake_atom[i][0];
|
|
buf[m++] = shake_atom[i][1];
|
|
buf[m++] = shake_atom[i][2];
|
|
buf[m++] = shake_atom[i][3];
|
|
buf[m++] = shake_type[i][0];
|
|
buf[m++] = shake_type[i][1];
|
|
buf[m++] = shake_type[i][2];
|
|
}
|
|
return m;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
unpack values in local atom-based arrays from exchange with another proc
|
|
------------------------------------------------------------------------- */
|
|
|
|
int FixShake::unpack_exchange(int nlocal, double *buf)
|
|
{
|
|
int m = 0;
|
|
int flag = shake_flag[nlocal] = static_cast<int> (buf[m++]);
|
|
if (flag == 1) {
|
|
shake_atom[nlocal][0] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][1] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][2] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][0] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][1] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][2] = static_cast<int> (buf[m++]);
|
|
} else if (flag == 2) {
|
|
shake_atom[nlocal][0] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][1] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][0] = static_cast<int> (buf[m++]);
|
|
} else if (flag == 3) {
|
|
shake_atom[nlocal][0] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][1] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][2] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][0] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][1] = static_cast<int> (buf[m++]);
|
|
} else if (flag == 4) {
|
|
shake_atom[nlocal][0] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][1] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][2] = static_cast<int> (buf[m++]);
|
|
shake_atom[nlocal][3] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][0] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][1] = static_cast<int> (buf[m++]);
|
|
shake_type[nlocal][2] = static_cast<int> (buf[m++]);
|
|
}
|
|
return m;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
enforce SHAKE constraints from rRESPA
|
|
prediction portion is different than Verlet
|
|
rRESPA updating of atom coords is done with full v, but only portions of f
|
|
------------------------------------------------------------------------- */
|
|
|
|
void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop)
|
|
{
|
|
// call stats only on outermost level
|
|
|
|
if (ilevel == nlevels_respa-1 && update->ntimestep == next_output) stats();
|
|
|
|
// perform SHAKE on every loop iteration of every rRESPA level
|
|
// except last loop iteration of inner levels
|
|
|
|
if (ilevel < nlevels_respa-1 && iloop == loop_respa[ilevel]-1) return;
|
|
|
|
// xshake = atom coords after next x update in innermost loop
|
|
// depends on rRESPA level
|
|
// for levels > 0 this includes more than one velocity update
|
|
// xshake = predicted position from call to this routine at level N =
|
|
// x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0)
|
|
|
|
double ***f_level = ((FixRespa *) modify->fix[ifix_respa])->f_level;
|
|
dtfsq = dtf_inner * step_respa[ilevel];
|
|
|
|
double invmass,dtfmsq;
|
|
int jlevel;
|
|
|
|
for (int i = 0; i < nlocal; i++) {
|
|
if (shake_flag[i]) {
|
|
invmass = 1.0 / mass[type[i]];
|
|
dtfmsq = dtfsq * invmass;
|
|
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
|
|
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
|
|
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
|
|
for (jlevel = 0; jlevel < ilevel; jlevel++) {
|
|
dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass;
|
|
xshake[i][0] += dtfmsq*f_level[i][jlevel][0];
|
|
xshake[i][1] += dtfmsq*f_level[i][jlevel][1];
|
|
xshake[i][2] += dtfmsq*f_level[i][jlevel][2];
|
|
}
|
|
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
|
|
}
|
|
|
|
// communicate results if necessary
|
|
|
|
if (nprocs > 1) comm->comm_fix(this);
|
|
|
|
// zero out SHAKE contribution to virial
|
|
|
|
vflag = vflag_in;
|
|
if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0;
|
|
|
|
// loop over clusters
|
|
|
|
int m;
|
|
for (int i = 0; i < nlist; i++) {
|
|
m = list[i];
|
|
if (shake_flag[m] == 2) shake2(m);
|
|
else if (shake_flag[m] == 3) shake3(m);
|
|
else if (shake_flag[m] == 4) shake4(m);
|
|
else shake3angle(m);
|
|
}
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
int FixShake::pack_comm(int n, int *list, double *buf, int *pbc_flags)
|
|
{
|
|
int i,j,m;
|
|
|
|
m = 0;
|
|
if (pbc_flags[0] == 0) {
|
|
for (i = 0; i < n; i++) {
|
|
j = list[i];
|
|
buf[m++] = xshake[j][0];
|
|
buf[m++] = xshake[j][1];
|
|
buf[m++] = xshake[j][2];
|
|
}
|
|
} else {
|
|
for (i = 0; i < n; i++) {
|
|
j = list[i];
|
|
buf[m++] = xshake[j][0] + pbc_flags[1]*domain->xprd;
|
|
buf[m++] = xshake[j][1] + pbc_flags[2]*domain->yprd;
|
|
buf[m++] = xshake[j][2] + pbc_flags[3]*domain->zprd;
|
|
}
|
|
}
|
|
return 3;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void FixShake::unpack_comm(int n, int first, double *buf)
|
|
{
|
|
int i,m,last;
|
|
|
|
m = 0;
|
|
last = first + n;
|
|
for (i = first; i < last; i++) {
|
|
xshake[i][0] = buf[m++];
|
|
xshake[i][1] = buf[m++];
|
|
xshake[i][2] = buf[m++];
|
|
}
|
|
}
|