diff --git a/src/USER-MISC/fix_filter_corotate.cpp b/src/USER-MISC/fix_filter_corotate.cpp
index 3e427f1f92533e858d3bcc6ef39c7870aad00cd5..7279859edb5fa301ced6705c230e3bb04ffde1d3 100644
--- a/src/USER-MISC/fix_filter_corotate.cpp
+++ b/src/USER-MISC/fix_filter_corotate.cpp
@@ -1,2065 +1,2062 @@
-/* ----------------------------------------------------------------------
- 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.
- ------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing author: Lukas Fath (KIT)
- some subroutines are from fix_shake.cpp
- ------------------------------------------------------------------------- */
-
-#include <mpi.h>
-#include "string.h"
-#include "stdlib.h"
-#include "fix_filter_corotate.h"
-#include "atom.h"
-#include "atom_vec.h"
-#include "molecule.h"
-#include "bond.h"
-#include "angle.h"
-#include "math_const.h"
-#include "update.h"
-#include "modify.h"
-#include "domain.h"
-#include "region.h"
-#include "memory.h"
-#include "error.h"
-#include "force.h"
-#include "comm.h"
-#include "error.h"
-#include "memory.h"
-#include "domain.h"
-#include "integrate.h"
-#include "respa.h"
-#include "timer.h"
-#include "neighbor.h"
-#include "citeme.h"
-
-using namespace LAMMPS_NS;
-using namespace MathConst;
-using namespace FixConst;
-
-// allocate space for static class variable
-
-FixFilterCorotate *FixFilterCorotate::fsptr;
-
-#define BIG 1.0e20
-#define MASSDELTA 0.1
-
-static const char cite_filter_corotate[] =
- "Mollified Impulse Method with Corotational Filter:\n\n"
- "@Article{Fath2017,\n"
- " Title = {A fast mollified impulse method for biomolecular atomistic simulations},\n"
- " Author = {L. Fath and M. Hochbruck and C.V. Singh},\n"
- " Journal = {Journal of Computational Physics},\n"
- " Year = {2017},\n"
- " Pages = {180 - 198},\n"
- " Volume = {333},\n\n"
- " Doi = {http://dx.doi.org/10.1016/j.jcp.2016.12.024},\n"
- " ISSN = {0021-9991},\n"
- " Keywords = {Mollified impulse method},\n"
- " Url = {http://www.sciencedirect.com/science/article/pii/S0021999116306787}\n"
- "}\n\n";
-
-/* ---------------------------------------------------------------------- */
-
-FixFilterCorotate::FixFilterCorotate(LAMMPS *lmp, int narg, char **arg):Fix(lmp,narg,arg)
-{
- if (lmp->citeme) lmp->citeme->add(cite_filter_corotate);
-
- MPI_Comm_rank(world,&me);
- MPI_Comm_size(world,&nprocs);
-
- molecular = atom->molecular;
- if (molecular == 0)
- error->all(FLERR,"Cannot use fix filter/corotate with non-molecular system");
- if (molecular == 2)
- error->all(FLERR,"Cannot use fix filter/corotate with molecule template system");
-
- // parse 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 = 3;
- 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(FLERR);
-
- // break if keyword that is not b,a,t,m
-
- } else if (isalpha(arg[next][0])) break;
-
- // read numeric args of b,a,t,m
-
- else if (mode == 'b') {
- int i = force->inumeric(FLERR,arg[next]);
- if (i < 1 || i > atom->nbondtypes)
- error->all(FLERR,"Invalid bond type index for fix filter/corotate");
- bond_flag[i] = 1;
-
- } else if (mode == 'a') {
- int i = force->inumeric(FLERR,arg[next]);
- if (i < 1 || i > atom->nangletypes)
- error->all(FLERR,"Invalid angle type index for fix filter/corotate");
- angle_flag[i] = 1;
-
- } else if (mode == 't') {
- int i = force->inumeric(FLERR,arg[next]);
- if (i < 1 || i > atom->ntypes)
- error->all(FLERR,"Invalid atom type index for fix filter/corotate");
- type_flag[i] = 1;
-
- } else if (mode == 'm') {
- double massone = force->numeric(FLERR,arg[next]);
- if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix filter/corotate");
- if (nmass == atom->ntypes)
- error->all(FLERR,"Too many masses for fix filter/corotate");
- mass_list[nmass++] = massone;
-
- } else error->all(FLERR,"Illegal fix filter/corotate command");
- next++;
- }
-
-
- // parse optional args
- nlevels_respa = 1; //initialize
-
- // 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];
-
- //grow_arrays
- array_atom = NULL;
- shake_flag = NULL;
- shake_atom = NULL;
- shake_type = NULL;
-
- grow_arrays(atom->nmax);
- atom->add_callback(0); //calls grow_arrays
-
- x_store = NULL;
-
-
- //STUFF
- g = NULL;
- help2 = NULL;
-
- dn1dx = dn2dx = dn3dx = NULL;
- n1 = n2 = n3 = del1 = del2 = del3 = NULL;
-
- memory->grow(help2,15,15,"FilterCorotate:help2");
- memory->grow(n1,3,"FilterCorotate:n1");
- memory->grow(n2,3,"FilterCorotate:n2");
- memory->grow(n3,3,"FilterCorotate:n3");
- memory->grow(del1,3,"FilterCorotate:del1");
- memory->grow(del2,3,"FilterCorotate:del2");
- memory->grow(del3,3,"FilterCorotate:del3");
- memory->grow(dn1dx,3,15,"FilterCorotate:dn1dx");
- memory->grow(dn2dx,3,15,"FilterCorotate:dn2dx");
- memory->grow(dn3dx,3,15,"FilterCorotate:dn3dx");
-
- memory->grow(g,15,"FilterCorotate:g");
-
- comm_forward = 3;
-
- // identify all clusters
-
- find_clusters();
-
- // initialize list of clusters to constrain
-
- maxlist = 0;
- list = NULL;
- clist_derv = NULL;
- clist_q0 = NULL; //list for derivative and ref. config
- clist_nselect1 = clist_nselect2 = NULL;
- clist_select1 = clist_select2 = NULL;
-
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixFilterCorotate::~FixFilterCorotate()
-{
- memory->destroy(g);
- memory->destroy(help2);
- memory->destroy(n1);
- memory->destroy(n2);
- memory->destroy(n3);
- memory->destroy(del1);
- memory->destroy(del2);
- memory->destroy(del3);
- memory->destroy(dn1dx);
- memory->destroy(dn2dx);
- memory->destroy(dn3dx);
-
- atom->delete_callback(id,2);
- atom->delete_callback(id,0);
-
- // delete locally stored arrays
-
- memory->destroy(shake_flag);
- memory->destroy(shake_atom);
- memory->destroy(shake_type);
- memory->destroy(array_atom);
-
- delete [] bond_flag;
- delete [] angle_flag;
- delete [] type_flag;
- delete [] mass_list;
-
- delete [] bond_distance;
- delete [] angle_distance;
-
- memory->destroy(list);
- memory->destroy(clist_derv);
- memory->destroy(clist_q0);
- memory->destroy(clist_nselect1);
- memory->destroy(clist_nselect2);
- memory->destroy(clist_select1);
- memory->destroy(clist_select2);
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-int FixFilterCorotate::setmask()
-{
- int mask = 0;
-
- mask |= PRE_NEIGHBOR;
- mask |= PRE_FORCE_RESPA;
- mask |= POST_FORCE_RESPA;
-
- return mask;
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-void FixFilterCorotate::init()
-{
- int i;
- // error if more than one filter
- int count = 0;
- for (i = 0; i < modify->nfix; i++)
- if (strcmp(modify->fix[i]->style,"filter/corotate") == 0) count++;
- if (count > 1) error->all(FLERR,"More than one fix filter/corotate");
-
- // check for fix shake:
- count = 0;
- for (i = 0; i < modify->nfix; i++)
- if (strcmp(modify->fix[i]->style,"shake") == 0) count++;
- if (count > 1) error->one(FLERR,"Fix shake and fix filter/corotate detected, are you sure?");
-
- // if rRESPA, find associated fix that must exist
- // could have changed locations in fix list since created
- // set ptrs to rRESPA variables
-
- if (strstr(update->integrate_style,"respa")) {
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
- }
- else error->all(FLERR,"Fix filter/corotate requires rRESPA!");
-
- // set equilibrium bond distances
-
- if (force->bond == NULL)
- error->all(FLERR,"Bond potential must be defined for fix filter/corotate");
- for (i = 1; i <= atom->nbondtypes; i++)
- bond_distance[i] = force->bond->equilibrium_distance(i);
-
- // set equilibrium angle value
-
- for (i = 1; i <= atom->nangletypes; i++) {
- angle_distance[i] = force->angle->equilibrium_angle(i);
- }
-}
-
-/* ----------------------------------------------------------------------
- * pre-integrator setup
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::setup_pre_neighbor()
-{
- pre_neighbor();
-}
-
-
-/* ----------------------------------------------------------------------
- * build list of clusters to constrain
- * if one or more atoms in cluster are on this proc,
- * this proc lists the cluster exactly once
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::pre_neighbor()
-{
- int atom1,atom2,atom3,atom4,atom5;
-
- // local copies of atom quantities
- // used until next re-neighboring
-
- x = atom->x;
- v = atom->v;
- f = atom->f;
- mass = atom->mass;
- rmass = atom->rmass;
- type = atom->type;
- nlocal = atom->nlocal;
-
- // extend size of list if necessary
-
- if (nlocal > maxlist) {
- maxlist = nlocal;
-
- memory->destroy(list);
- memory->destroy(clist_derv);
- memory->destroy(clist_q0);
- memory->destroy(clist_nselect1);
- memory->destroy(clist_nselect2);
- memory->destroy(clist_select1);
- memory->destroy(clist_select2);
-
- memory->create(list,maxlist,"FilterCorotate:list");
- memory->create(clist_derv,maxlist,15,15,"FilterCorotate:derivative");
- memory->create(clist_q0,maxlist,15,"FilterCorotate:q0");
- memory->create(clist_nselect1,maxlist,"FilterCorotate:nselect1");
- memory->create(clist_nselect2,maxlist,"FilterCorotate:nselect2");
- memory->create(clist_select1,maxlist,5,"FilterCorotate:select1");
- memory->create(clist_select2,maxlist,5,"FilterCorotate:select2");
- }
-
- // build list of 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,"Cluster atoms " TAGINT_FORMAT " " TAGINT_FORMAT
- " missing on proc %d at step " BIGINT_FORMAT,
- shake_atom[i][0],shake_atom[i][1],me,update->ntimestep);
- error->one(FLERR,str);
- }
- if (i <= atom1 && i <= atom2) list[nlist++] = i;
- } else if (shake_flag[i] == 1 || shake_flag[i] == 3) {
- 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,"Cluster atoms "
- TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
- " missing on proc %d at step " BIGINT_FORMAT,
- shake_atom[i][0],shake_atom[i][1],shake_atom[i][2],
- me,update->ntimestep);
- error->one(FLERR,str);
- }
- if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
- } else if (shake_flag[i] == 4) {
- 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,"Cluster atoms "
- TAGINT_FORMAT " " TAGINT_FORMAT " "
- TAGINT_FORMAT " " TAGINT_FORMAT
- " missing on proc %d at step " BIGINT_FORMAT,
- shake_atom[i][0],shake_atom[i][1],
- shake_atom[i][2],shake_atom[i][3],
- me,update->ntimestep);
- error->one(FLERR,str);
- }
- if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
- list[nlist++] = i;
- } else if (shake_flag[i] == 5) {
- 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]);
- atom5 = atom->map(shake_atom[i][4]);
- if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1 || atom5 == -1) {
- char str[128];
- sprintf(str,"Cluster atoms "
- TAGINT_FORMAT " " TAGINT_FORMAT " "
- TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
- " missing on proc %d at step " BIGINT_FORMAT,
- shake_atom[i][0],shake_atom[i][1],
- shake_atom[i][2],shake_atom[i][3],shake_atom[i][4],
- me,update->ntimestep);
- error->one(FLERR,str);
- }
- if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4 && i <= atom5)
- list[nlist++] = i;
- }
- }
-
- //also build reference config lists
- int m,N,oxy;
- double r0,r1,r2,r3,theta0;
- double m1,m2,m3,m4,m5,m_all;
-
- double *mass = atom->mass;
- int *type = atom->type;
-
- for (int i=0; i<nlist; i++)
- {
- m = list[i];
- N = shake_flag[m];
-
- //switch cluster type 3 to angle cluster:
- if (N == 3)
- {
- //make it an angle cluster:
- if (shake_type[m][2] == 0)
- shake_type[m][2] = angletype_findset(atom->map(shake_atom[m][0]),shake_atom[m][1],shake_atom[m][2],0);
- }
-
- if (N == 1) N = 3; //angle cluster
-
- if (N == 2) //cluster of size 2:
- {
- atom1 = atom->map(shake_atom[m][0]);
- atom2 = atom->map(shake_atom[m][1]);
-
- r0 = bond_distance[shake_type[m][0]];
-
- m1 = mass[type[atom1]];
- m2 = mass[type[atom2]];
- m_all = m1 + m2;
-
- double a1 = -m2/(m1+m2)*r0;
- double a2 = m1/(m1+m2)*r0;
-
- clist_q0[i][0] = a1;
- clist_q0[i][1] = 0;
- clist_q0[i][2] = 0;
-
- clist_q0[i][3] = a2;
- clist_q0[i][4] = 0;
- clist_q0[i][5] = 0;
-
- clist_nselect1[i] = 1;
- clist_select1[i][0] = 1;
-
- clist_nselect2[i] = 0;
- clist_select2[i][0] = 0;
- } else if (N == 3) //angle cluster
- {
- oxy = atom->map(shake_atom[m][0]);
- atom1 = atom->map(shake_atom[m][1]);
- atom2 = atom->map(shake_atom[m][2]);
-
- r0 = bond_distance[shake_type[m][0]];
- r1 = bond_distance[shake_type[m][1]];
-
-
- theta0 = angle_distance[shake_type[m][2]];
-
- m1 = mass[type[oxy]];
- m2 = mass[type[atom1]];
- m3 = mass[type[atom2]];
- m_all = m1 + m2 + m3;
-
- double alpha1 = 0.5*(MY_PI - theta0);
-
- double xcenter = m2/m_all*r0*sin(alpha1) + m3/m_all*r1*sin(alpha1);
- double ycenter = -m2/m_all*r0*cos(alpha1)+ m3/m_all*r1*cos(alpha1);
-
- double q1 = -xcenter;
- double q2 = -ycenter;
-
- double q4 = r0*sin(alpha1)-xcenter;
- double q5 = r0*cos(alpha1)-ycenter;
-
- double q7 = r1*sin(alpha1)-xcenter;
- double q8 = -r1*cos(alpha1)-ycenter;
-
- clist_q0[i][0] = q1;
- clist_q0[i][1] = q2;
- clist_q0[i][2] = 0;
-
- clist_q0[i][3] = q4;
- clist_q0[i][4] = q5;
- clist_q0[i][5] = 0;
-
- clist_q0[i][6] = q7;
- clist_q0[i][7] = q8;
- clist_q0[i][8] = 0;
-
- clist_nselect1[i] = 2;
- clist_select1[i][0] = 1; clist_select1[i][1] = 2;
- clist_nselect2[i] = 1;
- clist_select2[i][0] = 1;
- } else if (N == 4)
- {
- oxy = atom->map(shake_atom[m][0]);
- atom1 = atom->map(shake_atom[m][1]);
- atom2 = atom->map(shake_atom[m][2]);
- atom3 = atom->map(shake_atom[m][3]);
-
- r0 = bond_distance[shake_type[m][0]];
- r1 = bond_distance[shake_type[m][1]];
- r2 = bond_distance[shake_type[m][2]];
-
- m1 = atom->mass[atom->type[oxy]];
- m2 = atom->mass[atom->type[atom1]];
- m3 = atom->mass[atom->type[atom2]];
- m4 = atom->mass[atom->type[atom3]];
- m_all = m1 + m2 + m3 + m4;
-
- //how to get these angles?
- double alpha1 = MY_PI/2.57; //roughly 70 degrees
- double alpha2 = MY_PI/3;
- //ENSURE ycenter,zcenter = 0!
- //approximate xcenter, if r0 !=r1 != r2, exact if r0 = r1 = r2
- double xcenter = (m2*r0+m3*r1+m4*r2)/m_all*cos(alpha1);
-
- clist_q0[i][0] = -xcenter;
- clist_q0[i][1] = 0.0;
- clist_q0[i][2] = 0.0;
- clist_q0[i][3] = r0*cos(alpha1)-xcenter;
- clist_q0[i][4] = -r0*sin(alpha1);
- clist_q0[i][5] = 0.0;
- clist_q0[i][6] = r1*cos(alpha1)-xcenter;
- clist_q0[i][7] = r1*sin(alpha1)*cos(alpha2);
- clist_q0[i][8] = -r1*sin(alpha1)*sin(alpha2);
- clist_q0[i][9] = r2*cos(alpha1)-xcenter;
- clist_q0[i][10] = r2*sin(alpha1)*cos(alpha2);
- clist_q0[i][11] = r2*sin(alpha1)*sin(alpha2);
-
- clist_nselect1[i] = 3;
- clist_nselect2[i] = 2;
-
- clist_select1[i][0] = 1;clist_select1[i][1] = 2;clist_select1[i][2] = 3;
- clist_select2[i][0] = 2;clist_select2[i][1] = 3;
-
- //signum ensures correct ordering of three satellites
- //signum = sign(cross(x2-x1,x3-x1))T*(x1-x0))
- double del1[3], del2[3], del3[3];
- del1[0] = x[atom1][0]-x[oxy][0];
- del1[1] = x[atom1][1]-x[oxy][1];
- del1[2] = x[atom1][2]-x[oxy][2];
- domain->minimum_image(del1);
-
- del2[0] = x[atom2][0]-x[atom1][0];
- del2[1] = x[atom2][1]-x[atom1][1];
- del2[2] = x[atom2][2]-x[atom1][2];
- domain->minimum_image(del2);
-
- del3[0] = x[atom3][0]-x[atom1][0];
- del3[1] = x[atom3][1]-x[atom1][1];
- del3[2] = x[atom3][2]-x[atom1][2];
- domain->minimum_image(del3);
-
- double a = (del2[1])*(del3[2]) - (del2[2])*(del3[1]);
- double b = (del2[2])*(del3[0]) - (del2[0])*(del3[2]);
- double c = (del2[0])*(del3[1]) - (del2[1])*(del3[0]);
- int signum = sgn(a*(del1[0]) + b*(del1[1]) + c*(del1[2]));
-
- if (fabs(signum)!= 1) error->all(FLERR,"Wrong orientation in cluster of size 4 in fix filter/corotate!");
- clist_q0[i][8] *= signum;
- clist_q0[i][11] *= signum;
-
- } else if (N == 5)
- {
- oxy = atom->map(shake_atom[m][0]);
- atom1 = atom->map(shake_atom[m][1]);
- atom2 = atom->map(shake_atom[m][2]);
- atom3 = atom->map(shake_atom[m][3]);
- int c1 = atom->map(shake_atom[m][4]);
-
- r1 = bond_distance[shake_type[m][3]];
- r0 = bond_distance[shake_type[m][0]];
- theta0 = angle_distance[shake_type[m][2]];
-
- m1 = atom->mass[atom->type[oxy]];
- m2 = atom->mass[atom->type[atom1]];
- m3 = atom->mass[atom->type[atom2]];
- m4 = atom->mass[atom->type[atom3]];
- m5 = atom->mass[atom->type[c1]];
- m_all = m1 + m2 + m3 + m4 + m5;
-
- double alpha1 = MY_PI/2.57; //roughly 70 degrees
- double alpha2 = MY_PI/3;
- //ENSURE ycenter,zcenter = 0!
- double xcenter = -(m2+m3+m4)/m_all*r0*cos(alpha1) +r1*m5/m_all;
-
- clist_q0[i][0] = -xcenter;
- clist_q0[i][1] = 0.0;
- clist_q0[i][2] = 0.0;
- clist_q0[i][3] = -r0*cos(alpha1)-xcenter;
- clist_q0[i][4] = -r0*sin(alpha1);
- clist_q0[i][5] = 0.0;
- clist_q0[i][6] = -r0*cos(alpha1)-xcenter;
- clist_q0[i][7] = r0*sin(alpha1)*cos(alpha2);
- clist_q0[i][8] = r0*sin(alpha1)*sin(alpha2);
- clist_q0[i][9] = -r0*cos(alpha1)-xcenter;
- clist_q0[i][10] = r0*sin(alpha1)*cos(alpha2);
- clist_q0[i][11] = -r0*sin(alpha1)*sin(alpha2);
- clist_q0[i][12] = r1-xcenter;
- clist_q0[i][13] = 0.0;
- clist_q0[i][14] = 0.0;
- clist_nselect1[i] = 1;
- clist_nselect2[i] = 2;
-
- clist_select1[i][0] = 4;
- clist_select2[i][0] = 2;clist_select2[i][1] = 3;
-
- //signum ensures correct ordering of three satellites
- //signum = sign(cross(x2-x1,x3-x1))T*(x1-x0))
- double del1[3], del2[3], del3[3];
- del1[0] = x[atom1][0]-x[oxy][0];
- del1[1] = x[atom1][1]-x[oxy][1];
- del1[2] = x[atom1][2]-x[oxy][2];
- domain->minimum_image(del1);
-
- del2[0] = x[atom2][0]-x[atom1][0];
- del2[1] = x[atom2][1]-x[atom1][1];
- del2[2] = x[atom2][2]-x[atom1][2];
- domain->minimum_image(del2);
-
- del3[0] = x[atom3][0]-x[atom1][0];
- del3[1] = x[atom3][1]-x[atom1][1];
- del3[2] = x[atom3][2]-x[atom1][2];
- domain->minimum_image(del3);
-
- double a = (del2[1])*(del3[2]) - (del2[2])*(del3[1]);
- double b = (del2[2])*(del3[0]) - (del2[0])*(del3[2]);
- double c = (del2[0])*(del3[1]) - (del2[1])*(del3[0]);
- int signum = sgn(a*(del1[0]) + b*(del1[1]) + c*(del1[2]));
-
- if (fabs(signum)!= 1) error->all(FLERR,"Wrong orientation in cluster of size 5 in fix filter/corotate!");
- clist_q0[i][8] *= signum;
- clist_q0[i][11] *= signum;
- }
- else
- {
- error->all(FLERR,"Fix filter/corotate cluster with size > 5 not yet configured...");
- }
-
- }
-}
-
-/* ----------------------------------------------------------------------
- * setup, needs to patch missing setup_post_force_respa() in respa...
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::setup(int vflag)
-{
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
-}
-
-
-void FixFilterCorotate::setup_pre_force_respa(int vflag,int ilevel){
- pre_force_respa(vflag,ilevel,0);
-}
-
-//void FixFilterCorotate::setup_post_force_respa(int vflag,int ilevel){
-// post_force_respa(vflag,ilevel,0);
-//}
-
-
-double FixFilterCorotate::compute_array(int,int)
-{
- filter_inner();
- return 1;
-}
-
-void FixFilterCorotate::pre_force_respa(int vflag, int ilevel, int iloop)
-{
- if (ilevel == nlevels_respa-1)
- {
- filter_inner();
-
- //%Switch pointers x<-->x_filter
- x_store = atom->x;
- atom->x = array_atom;
- }
-}
-
-
-void FixFilterCorotate::post_force_respa(int vflag, int ilevel, int iloop)
-{
- if (ilevel == nlevels_respa-1)
- {
- atom->x = x_store;
-
- comm->forward_comm_fix(this); //comm forward of forces for outer filter
-
- filter_outer();
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixFilterCorotate::filter_inner()
-{
- //periodic boundaries, proc boundaries
-
- int nall = atom->nlocal + atom->nghost;
-
- double **x = atom->x;
-
- int nlocal = atom->nlocal;
- int* type = atom->type;
-
- //set array_atom to atom->x, in case some atoms are not part of any cluster:
- for (int i=0; i<nall; i++)
- {
- array_atom[i][0] = x[i][0];
- array_atom[i][1] = x[i][1];
- array_atom[i][2] = x[i][2];
- }
-
- int m;
- //iterate through list:
- for (int i=0; i<nlist; i++)
- {
- m = list[i];
-
- //filter:
- general_cluster(m,i);
- }
-
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-void FixFilterCorotate::filter_outer()
-{
- double sum1,sum2,sum3;
-
- double** f = atom->f;
- int nlocal = atom->nlocal;
-
- int m;
-
- for (int i=0; i<nlist; i++)
- {
- m = list[i];
- int dim = shake_flag[m];
- if (dim == 1) //case: angle
- dim = 3;
-
- for (int j = 0; j < dim; j++)
- {
- sum1 = sum2 = sum3 = 0;
- for (int k = 0; k < dim; k++)
- {
- //get local tag of k-th atom in cluster i:
- int kk = atom->map(shake_atom[m][k]);
- //apply derivative times force
- sum1 += clist_derv[i][3*j][3*k]*f[kk][0]+clist_derv[i][3*j][3*k+1]*f[kk][1]+clist_derv[i][3*j][3*k+2]*f[kk][2];
- sum2 += clist_derv[i][3*j+1][3*k]*f[kk][0]+clist_derv[i][3*j+1][3*k+1]*f[kk][1]+clist_derv[i][3*j+1][3*k+2]*f[kk][2];
- sum3 += clist_derv[i][3*j+2][3*k]*f[kk][0]+clist_derv[i][3*j+2][3*k+1]*f[kk][1]+clist_derv[i][3*j+2][3*k+2]*f[kk][2];
- }
-
- g[3*j] = sum1;
- g[3*j+1] = sum2;
- g[3*j+2] = sum3;
- }
-
- //replace force f with filtered value:
- for (int j = 0; j < dim; j++)
- {
- int kk = atom->map(shake_atom[m][j]);
-
- f[kk][0] = g[3*j];
- f[kk][1] = g[3*j+1];
- f[kk][2] = g[3*j+2];
- }
-
- }
-
-}
-
-
-/* ----------------------------------------------------------------------
- * 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 FixFilterCorotate::find_clusters()
-{
- int i,j,m,n,iatom;
- int flag,flag_all,nbuf,size;
- tagint tagprev;
- double massone;
- tagint *buf;
-
- if (me == 0 && screen)
- fprintf(screen,"Finding filter clusters ...\n");
-
-
-
- tagint *tag = atom->tag;
- int *type = atom->type;
- int *mask = atom->mask;
- double *mass = atom->mass;
- double *rmass = atom->rmass;
- int **nspecial = atom->nspecial;
- tagint **special = atom->special;
-
-
- int nlocal = atom->nlocal;
- int angles_allow = atom->avec->angles_allow;
-
- // 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_massflag[i][] = 1 if partner meets mass criterion, 0 if not
- // 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;
- memory->create(npartner,nlocal,"shake:npartner");
- memory->create(nshake,nlocal,"shake:nshake");
-
- tagint **partner_tag;
- int **partner_mask,**partner_type,**partner_massflag;
- int **partner_bondtype,**partner_shake,**partner_nshake;
- memory->create(partner_tag,nlocal,max,"shake:partner_tag");
- memory->create(partner_mask,nlocal,max,"shake:partner_mask");
- memory->create(partner_type,nlocal,max,"shake:partner_type");
- memory->create(partner_massflag,nlocal,max,"shake:partner_massflag");
- memory->create(partner_bondtype,nlocal,max,"shake:partner_bondtype");
- memory->create(partner_shake,nlocal,max,"shake:partner_shake");
- memory->create(partner_nshake,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_massflag, partner_bondtype
- // for bonded partners
- // requires communication for off-proc partners
- // -----------------------------------------------------
-
- // fill in mask, type, massflag, bondtype if own bond partner
- // info to store in buf for each off-proc bond = nper = 6
- // 2 atoms IDs in bond, space for mask, type, massflag, bondtype
- // nbufmax = largest buffer needed to hold info from any proc
-
- int nper = 6;
-
- 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_massflag[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];
- if (nmass) {
- if (rmass) massone = rmass[m];
- else massone = mass[type[m]];
- partner_massflag[i][j] = masscheck(massone);
- }
- n = bondtype_findset(i,tag[i],partner_tag[i][j],0);
- if (n) partner_bondtype[i][j] = n;
- else {
- n = bondtype_findset(m,tag[i],partner_tag[i][j],0);
- if (n) partner_bondtype[i][j] = n;
- }
- } else nbuf += nper;
- }
- }
-
- memory->create(buf,nbuf,"filter/corotate:buf");
-
- // 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;
- buf[size+4] = 0;
- n = bondtype_findset(i,tag[i],partner_tag[i][j],0);
- if (n) buf[size+5] = n;
- else buf[size+5] = 0;
- size += nper;
- }
- }
- }
-
- // cycle buffer around ring of procs back to self
-
- fsptr = this;
- comm->ring(size,sizeof(tagint),buf,1,ring_bonds,buf);
-
- // 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_massflag[i][j] = buf[m+4];
- partner_bondtype[i][j] = buf[m+5];
- m += nper;
- }
-
- memory->destroy(buf);
-
- // 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(FLERR,"Did not find fix filter/corotate 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
- // -----------------------------------------------------
-
- int np;
-
- for (i = 0; i < nlocal; i++) {
- nshake[i] = 0;
- np = npartner[i];
- for (j = 0; j < np; 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) {
- if (partner_massflag[i][j]) {
- partner_shake[i][j] = 1;
- nshake[i]++;
- continue;
- } else {
- if (rmass) massone = rmass[i];
- else massone = mass[type[i]];
- if (masscheck(massone)) {
- partner_shake[i][j] = 1;
- nshake[i]++;
- continue;
- }
- }
- }
- }
- }
-
- // -----------------------------------------------------
- // 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;
- }
- }
-
- memory->create(buf,nbuf,"filter/corotate:buf");
-
- // 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
-
- fsptr = this;
- comm->ring(size,sizeof(tagint),buf,2,ring_nshake,buf);
-
- // 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;
- }
-
- memory->destroy(buf);
-
- // -----------------------------------------------------
- // 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] > 4) flag = 1;
- MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
- if (flag_all) error->all(FLERR,"Cluster of more than 5 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(FLERR,"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_atom[i][4] = 0;
- shake_type[i][0] = 0;
- shake_type[i][1] = 0;
- shake_type[i][2] = 0;
- shake_type[i][3] = 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 && angles_allow) {
- n = angletype_findset(i,shake_atom[i][1],shake_atom[i][2],0);
- if (n <= 0) continue;
- //if (angle_flag[n]) { Take all angles!
- shake_flag[i] = 1;
- shake_type[i][2] = 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_atom[m][4] = shake_atom[i][4];
- shake_type[m][0] = shake_type[i][0];
- shake_type[m][1] = shake_type[i][1];
- shake_type[m][2] = shake_type[i][2];
- shake_type[m][3] = shake_type[i][3];
- } else nbuf += 11;
- }
- }
-
- memory->create(buf,nbuf,"filter/corotate:buf");
-
- // 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_atom[i][4];
- buf[size+7] = shake_type[i][0];
- buf[size+8] = shake_type[i][1];
- buf[size+9] = shake_type[i][2];
- buf[size+10] = shake_type[i][3];
- size += 11;
- }
- }
- }
-
- // cycle buffer around ring of procs back to self
-
- fsptr = this;
- comm->ring(size,sizeof(tagint),buf,3,ring_shake,NULL);
-
- memory->destroy(buf);
-
- // -----------------------------------------------------
- // free local memory
- // -----------------------------------------------------
-
- memory->destroy(npartner);
- memory->destroy(nshake);
- memory->destroy(partner_tag);
- memory->destroy(partner_mask);
- memory->destroy(partner_type);
- memory->destroy(partner_massflag);
- memory->destroy(partner_bondtype);
- memory->destroy(partner_shake);
- memory->destroy(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) {
- * bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
- * bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
- * angletype_findset(i,shake_atom[i][1],shake_atom[i][2],-1);
-} else if (shake_flag[i] == 2) {
- bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
-} else if (shake_flag[i] == 3) {
- bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
- bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
-} else if (shake_flag[i] == 4) {
- bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
- bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
- bondtype_findset(i,shake_atom[i][0],shake_atom[i][3],-1);
-}
-}*/
-
- // -----------------------------------------------------
- // print info on SHAKE clusters
- // -----------------------------------------------------
-
- int count1,count2,count3,count4,count5;
- count1 = count2 = count3 = count4 = count5 = 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++;
- else if (shake_flag[i] == 5) count5++;
- }
-
- 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);
- tmp = count5;
- MPI_Allreduce(&tmp,&count5,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 size 5 clusters\n",count5/5);
- 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 size 5 clusters\n",count5/5);
- fprintf(logfile," %d = # of frozen angles\n",count1/3);
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- * when receive buffer, scan bond partner IDs for atoms I own
- * if I own partner:
- * fill in mask and type and massflag
- * search for bond with 1st atom and fill in bondtype
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::ring_bonds(int ndatum, char *cbuf)
-{
- Atom *atom = fsptr->atom;
- double *rmass = atom->rmass;
- double *mass = atom->mass;
- int *mask = atom->mask;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int nmass = fsptr->nmass;
-
- tagint *buf = (tagint *) cbuf;
- int m,n;
- double massone;
-
- for (int i = 0; i < ndatum; i += 6) {
- m = atom->map(buf[i+1]);
- if (m >= 0 && m < nlocal) {
- buf[i+2] = mask[m];
- buf[i+3] = type[m];
- if (nmass) {
- if (rmass) massone = rmass[m];
- else massone = mass[type[m]];
- buf[i+4] = fsptr->masscheck(massone);
- }
- if (buf[i+5] == 0) {
- n = fsptr->bondtype_findset(m,buf[i],buf[i+1],0);
- if (n) buf[i+5] = n;
- }
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- * when receive buffer, scan bond partner IDs for atoms I own
- * if I own partner, fill in nshake value
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::ring_nshake(int ndatum, char *cbuf)
-{
- Atom *atom = fsptr->atom;
- int nlocal = atom->nlocal;
-
- int *nshake = fsptr->nshake;
-
- tagint *buf = (tagint *) cbuf;
- int m;
-
- for (int i = 0; i < ndatum; i += 3) {
- m = atom->map(buf[i+1]);
- if (m >= 0 && m < nlocal) buf[i+2] = nshake[m];
- }
-}
-
-/* ----------------------------------------------------------------------
- * when receive buffer, scan bond partner IDs for atoms I own
- * if I own partner, fill in nshake value
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::ring_shake(int ndatum, char *cbuf)
-{
- Atom *atom = fsptr->atom;
- int nlocal = atom->nlocal;
-
- int *shake_flag = fsptr->shake_flag;
- tagint **shake_atom = fsptr->shake_atom;
- int **shake_type = fsptr->shake_type;
-
- tagint *buf = (tagint *) cbuf;
- int m;
-
- for (int i = 0; i < ndatum; i += 11) {
- 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_atom[m][4] = buf[i+6];
- shake_type[m][0] = buf[i+7];
- shake_type[m][1] = buf[i+8];
- shake_type[m][2] = buf[i+9];
- shake_type[m][3] = buf[i+10];
- }
- }
-}
-
-/* ----------------------------------------------------------------------
- * check if massone is within MASSDELTA of any mass in mass_list
- * return 1 if yes, 0 if not
- * ------------------------------------------------------------------------- */
-
-int FixFilterCorotate::masscheck(double massone)
-{
- for (int i = 0; i < nmass; i++)
- if (fabs(mass_list[i]-massone) <= MASSDELTA) return 1;
- return 0;
-}
-
-/* --------------------------------------------------------------------------
- * filter one cluster of arbitrary size
- * ---------------------------------------------------------------------------*/
-
-void FixFilterCorotate::general_cluster(int index, int index_in_list)
-{
- //index: shake index, index_in_list: corresponding index in list
-
- //get q0, nselect:
- double *q0 = clist_q0[index_in_list];
- int nselect1 = clist_nselect1[index_in_list];
- int nselect2 = clist_nselect2[index_in_list];
- int *select1 = clist_select1[index_in_list];
- int *select2 = clist_select2[index_in_list];
-
- int N = shake_flag[index]; //number of atoms in cluster
- if (N == 1) N = 3; //angle cluster
-
- double**x = atom->x;
- double norm1, norm2, norm3;
-
- int* list_cluster = new int[N]; //contains local IDs of cluster atoms, 0 = center
- double* m = new double[N]; //contains local mass
- double *r = new double[N]; //contains r[i] = 1/||del[i]||
- double** del = new double*[N]; //contains del[i] = x_i-x_0
- for (int i = 0; i<N; i++)
- del[i] = new double[3];
-
- for (int i = 0; i < N; i++)
- {
- list_cluster[i] = atom->map(shake_atom[index][i]);
- m[i] = atom->mass[atom->type[list_cluster[i]]];
- }
-
- //%CALC r_i:
- for (int i = 1; i < N; i++)
- {
- del[i][0] = x[list_cluster[i]][0] - x[list_cluster[0]][0];
- del[i][1] = x[list_cluster[i]][1] - x[list_cluster[0]][1];
- del[i][2] = x[list_cluster[i]][2] - x[list_cluster[0]][2];
- domain->minimum_image(del[i]);
- r[i] = 1.0/sqrt(del[i][0]*del[i][0]+del[i][1]*del[i][1]+del[i][2]*del[i][2]);
- }
- //special case i=0: need del[0] for compact notation of array_atom later...
- del[0][0] = del[0][1] = del[0][2] = r[0] = 0.0;
-
- //calc n1:
- n1[0] = n1[1] = n1[2] = 0.0;
-
- int k;
- for (int i = 0; i < nselect1; i++)
- {
- k = select1[i];
- n1[0] += del[k][0]*r[k];
- n1[1] += del[k][1]*r[k];
- n1[2] += del[k][2]*r[k];
- }
- norm1 = 1.0/sqrt(n1[0]*n1[0]+n1[1]*n1[1]+n1[2]*n1[2]);
-
- n1[0] *= norm1;
- n1[1] *= norm1;
- n1[2] *= norm1;
-
- //calc n2:
- n2[0] = n2[1] = n2[2] = 0.0;
-
- for (int i = 0; i < nselect2; i++)
- {
- k = select2[i];
- n2[0] += del[k][0]*r[k];
- n2[1] += del[k][1]*r[k];
- n2[2] += del[k][2]*r[k];
- }
- if (!nselect2) //cluster of size 2, has only one direction
- n2[0] = n2[1] = n2[2] = 1.0;
-
- double alpha = n2[0]*n1[0] + n2[1]*n1[1] + n2[2]*n1[2];
-
- n2[0] -= alpha*n1[0];
- n2[1] -= alpha*n1[1];
- n2[2] -= alpha*n1[2];
-
- norm2 = 1.0/sqrt(n2[0]*n2[0]+n2[1]*n2[1]+n2[2]*n2[2]);
-
- n2[0] *= norm2;
- n2[1] *= norm2;
- n2[2] *= norm2;
-
- n3[0] = n1[1]*n2[2] - n1[2]*n2[1];
- n3[1] = n1[2]*n2[0] - n1[0]*n2[2];
- n3[2] = n1[0]*n2[1] - n1[1]*n2[0];
-
- norm3 = 1.0/sqrt(n3[0]*n3[0]+n3[1]*n3[1]+n3[2]*n3[2]);
-
- n3[0] *= norm3;
- n3[1] *= norm3;
- n3[2] *= norm3;
-
- //%x_filter:
- for (int i = 0; i < N; i++)
- {
- k = list_cluster[i];
- array_atom[k][0] = x[k][0]+q0[3*i]*n1[0]+q0[3*i+1]*n2[0]+q0[3*i+2]*n3[0];
- array_atom[k][1] = x[k][1]+q0[3*i]*n1[1]+q0[3*i+1]*n2[1]+q0[3*i+2]*n3[1];
- array_atom[k][2] = x[k][2]+q0[3*i]*n1[2]+q0[3*i+1]*n2[2]+q0[3*i+2]*n3[2];
- }
-
- double m_all = 0.0;
- for (int i = 0; i<N; i++)
- m_all += m[i];
-
- //%x+X_center
-
- for (int i = 0; i<N; i++)
- {
- k = list_cluster[i];
- for (int j = 0; j < N; j++)
- {
- array_atom[k][0] += m[j]/m_all*(del[j][0]-del[i][0]);
- array_atom[k][1] += m[j]/m_all*(del[j][1]-del[i][1]);
- array_atom[k][2] += m[j]/m_all*(del[j][2]-del[i][2]);
- }
- }
-
- //derivative:
- //dn1dx:
- double sum1[3][3*N];
- for (int i=0; i<3; i++)
- for (int j=0; j<3*N; j++)
- sum1[i][j] = 0;
-
- double I3mn1n1T[3][3]; //(I_3 - n1n1T)
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3; j++)
- I3mn1n1T[i][j] = -n1[i]*n1[j];
- I3mn1n1T[i][i] += 1.0;
- }
- // sum1 part of dn1dx:
- for (int l = 0; l < nselect1; l++)
- {
- k = select1[l];
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3; j++)
- {
- double help = del[k][i]*del[k][j]*r[k]*r[k]*r[k];
- sum1[i][j+3*k] -= help;
- sum1[i][j] += help;
- }
- sum1[i][i+3*k] += r[k];
- sum1[i][i] -= r[k];
- }
- }
- //dn1dx = norm1 * I3mn1n1T * sum1
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3*N; j++)
- {
- double sum = 0;
- for (int l = 0; l<3; l++)
- sum += I3mn1n1T[i][l]*sum1[l][j];
- dn1dx[i][j] = norm1*sum;
- }
- }
-
- //dn2dx: norm2 * I3mn2n2T * (I3mn1n1T*sum2 - rkn1pn1rk*dn1dx)
- double sum2[3][3*N];
- for (int i=0; i<3; i++)
- for (int j=0; j<3*N; j++)
- sum2[i][j] = 0;
-
- double I3mn2n2T[3][3]; //(I_3 - n2n2T)
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3; j++)
- I3mn2n2T[i][j] = -n2[i]*n2[j];
- I3mn2n2T[i][i] += 1.0;
- }
-
- // sum2 part of dn1dx:
- for (int l = 0; l < nselect2; l++)
- {
- k = select2[l];
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3; j++)
- {
- double help = del[k][i]*del[k][j]*r[k]*r[k]*r[k];
- sum2[i][j+3*k] -= help;
- sum2[i][j] += help;
- }
- sum2[i][i+3*k] += r[k];
- sum2[i][i] -= r[k];
- }
- }
-
- //prefaktor:
- double rkn1pn1rk[3][3];
- double rk[3]; rk[0] = rk[1] = rk[2] = 0.0;
-
- for (int i = 0; i < nselect2; i++)
- {
- k = select2[i];
- rk[0] += del[k][0]*r[k];
- rk[1] += del[k][1]*r[k];
- rk[2] += del[k][2]*r[k];
- }
- //rkn1pn1rk = rkT*n1*I3 + n1*rkT
- double scalar = rk[0]*n1[0]+rk[1]*n1[1]+rk[2]*n1[2];
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3; j++)
- rkn1pn1rk[i][j] = n1[i]*rk[j];
- rkn1pn1rk[i][i] += scalar;
- }
- //dn2dx: norm2 * I3mn2n2T * (I3mn1n1T*sum2 - rkn1pn1rk*dn1dx)
- //sum3 = (I3mn1n1T*sum2 - rkn1pn1rk*dn1dx)
- double sum3[3][3*N];
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3*N; j++)
- {
- double sum = 0;
- for (int l = 0; l<3; l++)
- sum += I3mn1n1T[i][l]*sum2[l][j] - rkn1pn1rk[i][l]*dn1dx[l][j];
- sum3[i][j] = sum;
- }
- }
- //dn2dx = norm2 * I3mn2n2T * sum3
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3*N; j++)
- {
- double sum = 0;
- for (int l = 0; l<3; l++)
- sum += I3mn2n2T[i][l]*sum3[l][j];
- dn2dx[i][j] = norm2*sum;
- }
- }
-
- //dn3dx = norm3 * I3mn3n3T * cross
- double I3mn3n3T[3][3]; //(I_3 - n3n3T)
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3; j++)
- I3mn3n3T[i][j] = -n3[i]*n3[j];
- I3mn3n3T[i][i] += 1.0;
- }
-
- double cross[3][3*N];
-
- for (int j=0; j<3*N; j++)
- {
- cross[0][j] = dn1dx[1][j]*n2[2] -dn1dx[2][j]*n2[1] + n1[1]*dn2dx[2][j]-n1[2]*dn2dx[1][j];
- cross[1][j] = dn1dx[2][j]*n2[0] -dn1dx[0][j]*n2[2] + n1[2]*dn2dx[0][j]-n1[0]*dn2dx[2][j];
- cross[2][j] = dn1dx[0][j]*n2[1] -dn1dx[1][j]*n2[0] + n1[0]*dn2dx[1][j]-n1[1]*dn2dx[0][j];
- }
-
- for (int i=0; i<3; i++)
- {
- for (int j=0; j<3*N; j++)
- {
- double sum = 0;
- for (int l = 0; l<3; l++)
- sum += I3mn3n3T[i][l]*cross[l][j];
- dn3dx[i][j] = norm3*sum;
- }
- }
- for (int l=0; l<N; l++)
- for (int i=0; i<3; i++)
- for (int j=0; j<3*N; j++)
- help2[i+3*l][j] = q0[3*l]*dn1dx[i][j] + q0[3*l+1]*dn2dx[i][j] + q0[3*l+2]*dn3dx[i][j];
-
- for (int j=0; j<N; j++)
- for (int l=0; l<N; l++)
- for (int i=0; i<3; i++)
- help2[i+3*l][i+3*j] += m[j]/m_all;
-
- //need transposed of help2:
- for (int ii = 0; ii < 3*N; ii++)
- for (int jj = 0; jj < 3*N; jj++)
- clist_derv[index_in_list][ii][jj] = help2[jj][ii];
-
- //free memory
- delete [] list_cluster;
- delete [] m;
- delete [] r;
- for (int i = 0; i<N; i++)
- delete [] del[i];
- delete [] del;
-}
-
-
-
-int FixFilterCorotate::pack_forward_comm(int n, int *list, double *buf,int pbc_flag, int *pbc)
-{
- int i,j,m;
- double**f = atom->f;
- m = 0;
- for (i = 0; i < n; i++)
- {
- j = list[i];
-
- buf[m++] = f[j][0];
- buf[m++] = f[j][1];
- buf[m++] = f[j][2];
- }
- return m;
-}
-
-void FixFilterCorotate::unpack_forward_comm(int n, int first, double *buf)
-{
- int i,m,last;
- double **f = atom->f;
- m = 0;
- last = first + n;
-
- for (i = first; i < last; i++)
- {
- f[i][0] = buf[m++];
- f[i][1] = buf[m++];
- f[i][2] = buf[m++];
- }
-}
-
-/* ----------------------------------------------------------------------
- * find a bond between global atom IDs n1 and n2 stored with local atom i
- * if find it:
- * if setflag = 0, return bond type
- * if setflag = -1/1, set bond type to negative/positive and return 0
- * if do not find it, return 0
- * ------------------------------------------------------------------------- */
-
-int FixFilterCorotate::bondtype_findset(int i, tagint n1, tagint n2, int setflag)
-{
- int m,nbonds;
- int *btype;
-
-
- tagint *tag = atom->tag;
- tagint **bond_atom = atom->bond_atom;
- nbonds = atom->num_bond[i];
-
- 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) {
- if (setflag == 0)
- return atom->bond_type[i][m];
-
- if ((setflag < 0 && atom->bond_type[i][m] > 0) ||
- (setflag > 0 && atom->bond_type[i][m] < 0))
- atom->bond_type[i][m] = -atom->bond_type[i][m];
-
- }
-
- return 0;
-}
-
-/* ----------------------------------------------------------------------
- * find an angle with global end atom IDs n1 and n2 stored with local atom i
- * if find it:
- * if setflag = 0, return angle type
- * if setflag = -1/1, set angle type to negative/positive and return 0
- * if do not find it, return 0
- * ------------------------------------------------------------------------- */
-
-int FixFilterCorotate::angletype_findset(int i, tagint n1, tagint n2, int setflag)
-{
- int m,nangles;
- int *atype;
-
-
- tagint **angle_atom1 = atom->angle_atom1;
- tagint **angle_atom3 = atom->angle_atom3;
- nangles = atom->num_angle[i];
-
- 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) {
- if (setflag == 0) {
- return atom->angle_type[i][m];
- }
-
- if ((setflag < 0 && atom->angle_type[i][m] > 0) ||
- (setflag > 0 && atom->angle_type[i][m] < 0))
- atom->angle_type[i][m] = -atom->angle_type[i][m];
-
- }
-
- return 0;
-}
-
-
-/* ----------------------------------------------------------------------
- * allocate local atom-based arrays
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::grow_arrays(int nmax)
-{
- memory->grow(array_atom,nmax,3,"FilterCorotate:peratomarray");
- memory->grow(shake_flag,nmax,"FilterCorotate::shake_flag");
- memory->grow(shake_atom,nmax,5,"FilterCorotate::shake_atom");
- memory->grow(shake_type,nmax,4,"FilterCorotate::shake_type");
-}
-
-/* ----------------------------------------------------------------------
- * memory usage of local atom-based arrays
- * ------------------------------------------------------------------------- */
-
-double FixFilterCorotate::memory_usage()
-{
- double bytes = 0;
- //GROW:
- bytes += 3*sizeof(double) + 5*sizeof(tagint) + 5*sizeof(int);
- //clist
- bytes += 13*atom->nlocal*sizeof(int);
- bytes += 15*16*nlocal*sizeof(double);
-
- //fixed:
- int nb = atom->nbondtypes+1;
- int na = atom->nangletypes+1;
- int nt = atom->ntypes+1;
- bytes += (nb+na+nt)*sizeof(int);
- bytes += (nt-1+nb+na+15*15+18+10*15)*sizeof(double);
-
- //output:
- //bytes += (2*nb+2*na)*sizeof(int);
- //bytes += (6*na+6*nb)*sizeof(double);
-
- return bytes;
-}
-
-
-/* ----------------------------------------------------------------------
- * copy values within local atom-based arrays
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::copy_arrays(int i, int j, int delflag)
-{
- 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];
- shake_type[j][2] = shake_type[i][2];
- } 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];
- } else if (flag == 5) {
- 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_atom[j][4] = shake_atom[i][4];
- shake_type[j][0] = shake_type[i][0];
- shake_type[j][1] = shake_type[i][1];
- shake_type[j][2] = shake_type[i][2];
- shake_type[j][3] = shake_type[i][3];
- }
-}
-
-/* ----------------------------------------------------------------------
- * initialize one atom's array values, called when atom is created
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::set_arrays(int i)
-{
- shake_flag[i] = 0;
-}
-
-/* ----------------------------------------------------------------------
- * update one atom's array values
- * called when molecule is created from fix gcmc
- * ------------------------------------------------------------------------- */
-
-void FixFilterCorotate::update_arrays(int i, int atom_offset)
-{
- int flag = shake_flag[i];
-
- if (flag == 1) {
- shake_atom[i][0] += atom_offset;
- shake_atom[i][1] += atom_offset;
- shake_atom[i][2] += atom_offset;
- } else if (flag == 2) {
- shake_atom[i][0] += atom_offset;
- shake_atom[i][1] += atom_offset;
- } else if (flag == 3) {
- shake_atom[i][0] += atom_offset;
- shake_atom[i][1] += atom_offset;
- shake_atom[i][2] += atom_offset;
- } else if (flag == 4) {
- shake_atom[i][0] += atom_offset;
- shake_atom[i][1] += atom_offset;
- shake_atom[i][2] += atom_offset;
- shake_atom[i][3] += atom_offset;
- } else if (flag == 5) {
- shake_atom[i][0] += atom_offset;
- shake_atom[i][1] += atom_offset;
- shake_atom[i][2] += atom_offset;
- shake_atom[i][3] += atom_offset;
- shake_atom[i][4] += atom_offset;
- }
-}
-
-
-/* ----------------------------------------------------------------------
- * pack values in local atom-based arrays for exchange with another proc
- * ------------------------------------------------------------------------- */
-
-int FixFilterCorotate::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];
- buf[m++] = shake_type[i][2];
- } 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];
- } else if (flag == 5) {
- 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_atom[i][4];
- buf[m++] = shake_type[i][0];
- buf[m++] = shake_type[i][1];
- buf[m++] = shake_type[i][2];
- buf[m++] = shake_type[i][3];
- }
- return m;
-}
-
-/* ----------------------------------------------------------------------
- * unpack values in local atom-based arrays from exchange with another proc
- * ------------------------------------------------------------------------- */
-
-int FixFilterCorotate::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<tagint> (buf[m++]);
- shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][2] = static_cast<tagint> (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<tagint> (buf[m++]);
- shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
- shake_type[nlocal][0] = static_cast<int> (buf[m++]);
- } else if (flag == 3) {
- shake_atom[nlocal][0] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][2] = static_cast<tagint> (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 == 4) {
- shake_atom[nlocal][0] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][2] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][3] = static_cast<tagint> (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 == 5) {
- shake_atom[nlocal][0] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][2] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][3] = static_cast<tagint> (buf[m++]);
- shake_atom[nlocal][4] = static_cast<tagint> (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++]);
- shake_type[nlocal][3] = static_cast<int> (buf[m++]);
- }
- return m;
-}
+/* ----------------------------------------------------------------------
+ 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.
+ ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Lukas Fath (KIT)
+ some subroutines are from fix_shake.cpp
+ ------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include <string.h>
+#include <stdlib.h>
+#include "fix_filter_corotate.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "molecule.h"
+#include "bond.h"
+#include "angle.h"
+#include "math_const.h"
+#include "update.h"
+#include "modify.h"
+#include "domain.h"
+#include "region.h"
+#include "memory.h"
+#include "error.h"
+#include "force.h"
+#include "comm.h"
+#include "error.h"
+#include "memory.h"
+#include "domain.h"
+#include "integrate.h"
+#include "respa.h"
+#include "neighbor.h"
+#include "citeme.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+using namespace FixConst;
+
+// allocate space for static class variable
+
+FixFilterCorotate *FixFilterCorotate::fsptr;
+
+#define BIG 1.0e20
+#define MASSDELTA 0.1
+
+static const char cite_filter_corotate[] =
+"Mollified Impulse Method with Corotational Filter:\n\n"
+"@Article{Fath2017,\n"
+" Title ="
+"{A fast mollified impulse method for biomolecular atomistic simulations},\n"
+" Author = {L. Fath and M. Hochbruck and C.V. Singh},\n"
+" Journal = {Journal of Computational Physics},\n"
+" Year = {2017},\n"
+" Pages = {180 - 198},\n"
+" Volume = {333},\n\n"
+" Doi = {http://dx.doi.org/10.1016/j.jcp.2016.12.024},\n"
+" ISSN = {0021-9991},\n"
+" Keywords = {Mollified impulse method},\n"
+" Url = {http://www.sciencedirect.com/science/article/pii/S0021999116306787}\n"
+"}\n\n";
+
+/* ---------------------------------------------------------------------- */
+
+FixFilterCorotate::FixFilterCorotate(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp,narg,arg)
+{
+ if (lmp->citeme) lmp->citeme->add(cite_filter_corotate);
+
+ MPI_Comm_rank(world,&me);
+ MPI_Comm_size(world,&nprocs);
+
+ molecular = atom->molecular;
+ if (molecular == 0)
+ error->all(FLERR,"Cannot use fix filter/corotate "
+ "with non-molecular system");
+ if (molecular == 2)
+ error->all(FLERR,"Cannot use fix filter/corotate "
+ "with molecular template system");
+
+ // parse 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 = 3;
+ 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(FLERR);
+
+ // break if keyword that is not b,a,t,m
+
+ } else if (isalpha(arg[next][0])) break;
+
+ // read numeric args of b,a,t,m
+
+ else if (mode == 'b') {
+ int i = force->inumeric(FLERR,arg[next]);
+ if (i < 1 || i > atom->nbondtypes)
+ error->all(FLERR,"Invalid bond type index for fix filter/corotate");
+ bond_flag[i] = 1;
+
+ } else if (mode == 'a') {
+ int i = force->inumeric(FLERR,arg[next]);
+ if (i < 1 || i > atom->nangletypes)
+ error->all(FLERR,"Invalid angle type index for fix filter/corotate");
+ angle_flag[i] = 1;
+
+ } else if (mode == 't') {
+ int i = force->inumeric(FLERR,arg[next]);
+ if (i < 1 || i > atom->ntypes)
+ error->all(FLERR,"Invalid atom type index for fix filter/corotate");
+ type_flag[i] = 1;
+
+ } else if (mode == 'm') {
+ double massone = force->numeric(FLERR,arg[next]);
+ if (massone == 0.0)
+ error->all(FLERR,"Invalid atom mass for fix filter/corotate");
+ if (nmass == atom->ntypes)
+ error->all(FLERR,"Too many masses for fix filter/corotate");
+ mass_list[nmass++] = massone;
+
+ } else error->all(FLERR,"Illegal fix filter/corotate command");
+ next++;
+ }
+
+
+ // parse optional args
+ nlevels_respa = 1; //initialize
+
+ // 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];
+
+ //grow_arrays
+ array_atom = NULL;
+ shake_flag = NULL;
+ shake_atom = NULL;
+ shake_type = NULL;
+
+ grow_arrays(atom->nmax);
+ atom->add_callback(0); //calls grow_arrays
+
+ x_store = NULL;
+
+ //STUFF
+ g = NULL;
+ help2 = NULL;
+
+ dn1dx = dn2dx = dn3dx = NULL;
+ n1 = n2 = n3 = del1 = del2 = del3 = NULL;
+
+ memory->grow(help2,15,15,"FilterCorotate:help2");
+ memory->grow(n1,3,"FilterCorotate:n1");
+ memory->grow(n2,3,"FilterCorotate:n2");
+ memory->grow(n3,3,"FilterCorotate:n3");
+ memory->grow(del1,3,"FilterCorotate:del1");
+ memory->grow(del2,3,"FilterCorotate:del2");
+ memory->grow(del3,3,"FilterCorotate:del3");
+ memory->grow(dn1dx,3,15,"FilterCorotate:dn1dx");
+ memory->grow(dn2dx,3,15,"FilterCorotate:dn2dx");
+ memory->grow(dn3dx,3,15,"FilterCorotate:dn3dx");
+
+ memory->grow(g,15,"FilterCorotate:g");
+
+ comm_forward = 3;
+
+ // identify all clusters
+
+ find_clusters();
+
+ // initialize list of clusters to constrain
+
+ maxlist = 0;
+ list = NULL;
+ clist_derv = NULL;
+ clist_q0 = NULL; //list for derivative and ref. config
+ clist_nselect1 = clist_nselect2 = NULL;
+ clist_select1 = clist_select2 = NULL;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixFilterCorotate::~FixFilterCorotate()
+{
+ memory->destroy(g);
+ memory->destroy(help2);
+ memory->destroy(n1);
+ memory->destroy(n2);
+ memory->destroy(n3);
+ memory->destroy(del1);
+ memory->destroy(del2);
+ memory->destroy(del3);
+ memory->destroy(dn1dx);
+ memory->destroy(dn2dx);
+ memory->destroy(dn3dx);
+
+ atom->delete_callback(id,2);
+ atom->delete_callback(id,0);
+
+ // delete locally stored arrays
+
+ memory->destroy(shake_flag);
+ memory->destroy(shake_atom);
+ memory->destroy(shake_type);
+ memory->destroy(array_atom);
+
+ delete [] bond_flag;
+ delete [] angle_flag;
+ delete [] type_flag;
+ delete [] mass_list;
+
+ delete [] bond_distance;
+ delete [] angle_distance;
+
+ memory->destroy(list);
+ memory->destroy(clist_derv);
+ memory->destroy(clist_q0);
+ memory->destroy(clist_nselect1);
+ memory->destroy(clist_nselect2);
+ memory->destroy(clist_select1);
+ memory->destroy(clist_select2);
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int FixFilterCorotate::setmask()
+{
+ int mask = 0;
+
+ mask |= PRE_NEIGHBOR;
+ mask |= PRE_FORCE_RESPA;
+ mask |= POST_FORCE_RESPA;
+
+ return mask;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixFilterCorotate::init()
+{
+ int i;
+ // error if more than one filter
+ int count = 0;
+ for (i = 0; i < modify->nfix; i++){
+ if (strcmp(modify->fix[i]->style,"filter/corotate") == 0) count++;
+ }
+ if (count > 1) error->all(FLERR,"More than one fix filter/corotate");
+
+ // check for fix shake:
+ count = 0;
+ for (i = 0; i < modify->nfix; i++){
+ if (strcmp(modify->fix[i]->style,"shake") == 0) count++;
+ }
+ if (count > 1)
+ error->one(FLERR,"Both fix shake and fix filter/corotate detected.");
+
+ // if rRESPA, find associated fix that must exist
+ // could have changed locations in fix list since created
+ // set ptrs to rRESPA variables
+
+ if (strstr(update->integrate_style,"respa")) {
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ }
+ else error->all(FLERR,"Fix filter/corotate requires rRESPA!");
+
+ // set equilibrium bond distances
+
+ if (force->bond == NULL)
+ error->all(FLERR,"Bond potential must be defined for fix filter/corotate");
+ for (i = 1; i <= atom->nbondtypes; i++)
+ bond_distance[i] = force->bond->equilibrium_distance(i);
+
+ // set equilibrium angle value
+
+ for (i = 1; i <= atom->nangletypes; i++) {
+ angle_distance[i] = force->angle->equilibrium_angle(i);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * pre-integrator setup
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::setup_pre_neighbor()
+{
+ pre_neighbor();
+}
+
+
+/* ----------------------------------------------------------------------
+ * build list of clusters to constrain
+ * if one or more atoms in cluster are on this proc,
+ * this proc lists the cluster exactly once
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::pre_neighbor()
+{
+ int atom1,atom2,atom3,atom4,atom5;
+
+ // local copies of atom quantities
+ // used until next re-neighboring
+
+ x = atom->x;
+ v = atom->v;
+ f = atom->f;
+ mass = atom->mass;
+ rmass = atom->rmass;
+ type = atom->type;
+ nlocal = atom->nlocal;
+
+ // extend size of list if necessary
+
+ if (nlocal > maxlist) {
+ maxlist = nlocal;
+
+ memory->destroy(list);
+ memory->destroy(clist_derv);
+ memory->destroy(clist_q0);
+ memory->destroy(clist_nselect1);
+ memory->destroy(clist_nselect2);
+ memory->destroy(clist_select1);
+ memory->destroy(clist_select2);
+
+ memory->create(list,maxlist,"FilterCorotate:list");
+ memory->create(clist_derv,maxlist,15,15,"FilterCorotate:derivative");
+ memory->create(clist_q0,maxlist,15,"FilterCorotate:q0");
+ memory->create(clist_nselect1,maxlist,"FilterCorotate:nselect1");
+ memory->create(clist_nselect2,maxlist,"FilterCorotate:nselect2");
+ memory->create(clist_select1,maxlist,5,"FilterCorotate:select1");
+ memory->create(clist_select2,maxlist,5,"FilterCorotate:select2");
+ }
+
+ // build list of 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,"Cluster atoms " TAGINT_FORMAT " " TAGINT_FORMAT
+ " missing on proc %d at step " BIGINT_FORMAT,
+ shake_atom[i][0],shake_atom[i][1],me,update->ntimestep);
+ error->one(FLERR,str);
+ }
+ if (i <= atom1 && i <= atom2) list[nlist++] = i;
+ } else if (shake_flag[i] == 1 || shake_flag[i] == 3) {
+ 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,"Cluster atoms "
+ TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
+ " missing on proc %d at step " BIGINT_FORMAT,
+ shake_atom[i][0],shake_atom[i][1],shake_atom[i][2],
+ me,update->ntimestep);
+ error->one(FLERR,str);
+ }
+ if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
+ } else if (shake_flag[i] == 4) {
+ 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,"Cluster atoms "
+ TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT " " TAGINT_FORMAT
+ " missing on proc %d at step " BIGINT_FORMAT,
+ shake_atom[i][0],shake_atom[i][1],
+ shake_atom[i][2],shake_atom[i][3],
+ me,update->ntimestep);
+ error->one(FLERR,str);
+ }
+ if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
+ list[nlist++] = i;
+ } else if (shake_flag[i] == 5) {
+ 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]);
+ atom5 = atom->map(shake_atom[i][4]);
+ if (atom1 == -1 || atom2 == -1 || atom3 == -1 ||
+ atom4 == -1 || atom5 == -1) {
+ char str[128];
+ sprintf(str,"Cluster atoms "
+ TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
+ " missing on proc %d at step " BIGINT_FORMAT,
+ shake_atom[i][0],shake_atom[i][1],
+ shake_atom[i][2],shake_atom[i][3],shake_atom[i][4],
+ me,update->ntimestep);
+ error->one(FLERR,str);
+ }
+ if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4 && i <= atom5)
+ list[nlist++] = i;
+ }
+ }
+ }
+
+ //also build reference config lists
+ int m,N,oxy;
+ double r0,r1,r2,theta0;
+ double m1,m2,m3,m4,m5,m_all;
+
+ double *mass = atom->mass;
+ int *type = atom->type;
+
+ for (int i=0; i<nlist; i++)
+ {
+ m = list[i];
+ N = shake_flag[m];
+
+ //switch cluster type 3 to angle cluster:
+ if (N == 3)
+ {
+ //make it an angle cluster:
+ if (shake_type[m][2] == 0)
+ shake_type[m][2] = angletype_findset(atom->map(shake_atom[m][0]),
+ shake_atom[m][1],shake_atom[m][2],0);
+ }
+
+ if (N == 1) N = 3; //angle cluster
+
+ if (N == 2) //cluster of size 2:
+ {
+ atom1 = atom->map(shake_atom[m][0]);
+ atom2 = atom->map(shake_atom[m][1]);
+
+ r0 = bond_distance[shake_type[m][0]];
+
+ m1 = mass[type[atom1]];
+ m2 = mass[type[atom2]];
+ m_all = m1 + m2;
+
+ double a1 = -m2/(m1+m2)*r0;
+ double a2 = m1/(m1+m2)*r0;
+
+ clist_q0[i][0] = a1;
+ clist_q0[i][1] = 0;
+ clist_q0[i][2] = 0;
+
+ clist_q0[i][3] = a2;
+ clist_q0[i][4] = 0;
+ clist_q0[i][5] = 0;
+
+ clist_nselect1[i] = 1;
+ clist_select1[i][0] = 1;
+
+ clist_nselect2[i] = 0;
+ clist_select2[i][0] = 0;
+ } else if (N == 3) //angle cluster
+ {
+ oxy = atom->map(shake_atom[m][0]);
+ atom1 = atom->map(shake_atom[m][1]);
+ atom2 = atom->map(shake_atom[m][2]);
+
+ r0 = bond_distance[shake_type[m][0]];
+ r1 = bond_distance[shake_type[m][1]];
+
+
+ theta0 = angle_distance[shake_type[m][2]];
+
+ m1 = mass[type[oxy]];
+ m2 = mass[type[atom1]];
+ m3 = mass[type[atom2]];
+ m_all = m1 + m2 + m3;
+
+ double alpha1 = 0.5*(MY_PI - theta0);
+
+ double xcenter = m2/m_all*r0*sin(alpha1) + m3/m_all*r1*sin(alpha1);
+ double ycenter = -m2/m_all*r0*cos(alpha1)+ m3/m_all*r1*cos(alpha1);
+
+ double q1 = -xcenter;
+ double q2 = -ycenter;
+
+ double q4 = r0*sin(alpha1)-xcenter;
+ double q5 = r0*cos(alpha1)-ycenter;
+
+ double q7 = r1*sin(alpha1)-xcenter;
+ double q8 = -r1*cos(alpha1)-ycenter;
+
+ clist_q0[i][0] = q1;
+ clist_q0[i][1] = q2;
+ clist_q0[i][2] = 0;
+
+ clist_q0[i][3] = q4;
+ clist_q0[i][4] = q5;
+ clist_q0[i][5] = 0;
+
+ clist_q0[i][6] = q7;
+ clist_q0[i][7] = q8;
+ clist_q0[i][8] = 0;
+
+ clist_nselect1[i] = 2;
+ clist_select1[i][0] = 1; clist_select1[i][1] = 2;
+ clist_nselect2[i] = 1;
+ clist_select2[i][0] = 1;
+ } else if (N == 4)
+ {
+ oxy = atom->map(shake_atom[m][0]);
+ atom1 = atom->map(shake_atom[m][1]);
+ atom2 = atom->map(shake_atom[m][2]);
+ atom3 = atom->map(shake_atom[m][3]);
+
+ r0 = bond_distance[shake_type[m][0]];
+ r1 = bond_distance[shake_type[m][1]];
+ r2 = bond_distance[shake_type[m][2]];
+
+ m1 = atom->mass[atom->type[oxy]];
+ m2 = atom->mass[atom->type[atom1]];
+ m3 = atom->mass[atom->type[atom2]];
+ m4 = atom->mass[atom->type[atom3]];
+ m_all = m1 + m2 + m3 + m4;
+
+ //how to get these angles?
+ double alpha1 = MY_PI/2.57; //roughly 70 degrees
+ double alpha2 = MY_PI/3;
+ //ENSURE ycenter,zcenter = 0!
+ //approximate xcenter, if r0 !=r1 != r2, exact if r0 = r1 = r2
+ double xcenter = (m2*r0+m3*r1+m4*r2)/m_all*cos(alpha1);
+
+ clist_q0[i][0] = -xcenter;
+ clist_q0[i][1] = 0.0;
+ clist_q0[i][2] = 0.0;
+ clist_q0[i][3] = r0*cos(alpha1)-xcenter;
+ clist_q0[i][4] = -r0*sin(alpha1);
+ clist_q0[i][5] = 0.0;
+ clist_q0[i][6] = r1*cos(alpha1)-xcenter;
+ clist_q0[i][7] = r1*sin(alpha1)*cos(alpha2);
+ clist_q0[i][8] = -r1*sin(alpha1)*sin(alpha2);
+ clist_q0[i][9] = r2*cos(alpha1)-xcenter;
+ clist_q0[i][10] = r2*sin(alpha1)*cos(alpha2);
+ clist_q0[i][11] = r2*sin(alpha1)*sin(alpha2);
+
+ clist_nselect1[i] = 3;
+ clist_nselect2[i] = 2;
+
+ clist_select1[i][0] = 1;clist_select1[i][1] = 2;clist_select1[i][2] = 3;
+ clist_select2[i][0] = 2;clist_select2[i][1] = 3;
+
+ //signum ensures correct ordering of three satellites
+ //signum = sign(cross(x2-x1,x3-x1))T*(x1-x0))
+ double del1[3], del2[3], del3[3];
+ del1[0] = x[atom1][0]-x[oxy][0];
+ del1[1] = x[atom1][1]-x[oxy][1];
+ del1[2] = x[atom1][2]-x[oxy][2];
+ domain->minimum_image(del1);
+
+ del2[0] = x[atom2][0]-x[atom1][0];
+ del2[1] = x[atom2][1]-x[atom1][1];
+ del2[2] = x[atom2][2]-x[atom1][2];
+ domain->minimum_image(del2);
+
+ del3[0] = x[atom3][0]-x[atom1][0];
+ del3[1] = x[atom3][1]-x[atom1][1];
+ del3[2] = x[atom3][2]-x[atom1][2];
+ domain->minimum_image(del3);
+
+ double a = (del2[1])*(del3[2]) - (del2[2])*(del3[1]);
+ double b = (del2[2])*(del3[0]) - (del2[0])*(del3[2]);
+ double c = (del2[0])*(del3[1]) - (del2[1])*(del3[0]);
+ int signum = sgn(a*(del1[0]) + b*(del1[1]) + c*(del1[2]));
+
+ if (fabs(signum)!= 1)
+ error->all(FLERR,"Wrong orientation in cluster of size 4"
+ "in fix filter/corotate!");
+ clist_q0[i][8] *= signum;
+ clist_q0[i][11] *= signum;
+
+ } else if (N == 5)
+ {
+ oxy = atom->map(shake_atom[m][0]);
+ atom1 = atom->map(shake_atom[m][1]);
+ atom2 = atom->map(shake_atom[m][2]);
+ atom3 = atom->map(shake_atom[m][3]);
+ int c1 = atom->map(shake_atom[m][4]);
+
+ r1 = bond_distance[shake_type[m][3]];
+ r0 = bond_distance[shake_type[m][0]];
+ theta0 = angle_distance[shake_type[m][2]];
+
+ m1 = atom->mass[atom->type[oxy]];
+ m2 = atom->mass[atom->type[atom1]];
+ m3 = atom->mass[atom->type[atom2]];
+ m4 = atom->mass[atom->type[atom3]];
+ m5 = atom->mass[atom->type[c1]];
+ m_all = m1 + m2 + m3 + m4 + m5;
+
+ double alpha1 = MY_PI/2.57; //roughly 70 degrees
+ double alpha2 = MY_PI/3;
+ //ENSURE ycenter,zcenter = 0!
+ double xcenter = -(m2+m3+m4)/m_all*r0*cos(alpha1) +r1*m5/m_all;
+
+ clist_q0[i][0] = -xcenter;
+ clist_q0[i][1] = 0.0;
+ clist_q0[i][2] = 0.0;
+ clist_q0[i][3] = -r0*cos(alpha1)-xcenter;
+ clist_q0[i][4] = -r0*sin(alpha1);
+ clist_q0[i][5] = 0.0;
+ clist_q0[i][6] = -r0*cos(alpha1)-xcenter;
+ clist_q0[i][7] = r0*sin(alpha1)*cos(alpha2);
+ clist_q0[i][8] = r0*sin(alpha1)*sin(alpha2);
+ clist_q0[i][9] = -r0*cos(alpha1)-xcenter;
+ clist_q0[i][10] = r0*sin(alpha1)*cos(alpha2);
+ clist_q0[i][11] = -r0*sin(alpha1)*sin(alpha2);
+ clist_q0[i][12] = r1-xcenter;
+ clist_q0[i][13] = 0.0;
+ clist_q0[i][14] = 0.0;
+ clist_nselect1[i] = 1;
+ clist_nselect2[i] = 2;
+
+ clist_select1[i][0] = 4;
+ clist_select2[i][0] = 2;clist_select2[i][1] = 3;
+
+ //signum ensures correct ordering of three satellites
+ //signum = sign(cross(x2-x1,x3-x1))T*(x1-x0))
+ double del1[3], del2[3], del3[3];
+ del1[0] = x[atom1][0]-x[oxy][0];
+ del1[1] = x[atom1][1]-x[oxy][1];
+ del1[2] = x[atom1][2]-x[oxy][2];
+ domain->minimum_image(del1);
+
+ del2[0] = x[atom2][0]-x[atom1][0];
+ del2[1] = x[atom2][1]-x[atom1][1];
+ del2[2] = x[atom2][2]-x[atom1][2];
+ domain->minimum_image(del2);
+
+ del3[0] = x[atom3][0]-x[atom1][0];
+ del3[1] = x[atom3][1]-x[atom1][1];
+ del3[2] = x[atom3][2]-x[atom1][2];
+ domain->minimum_image(del3);
+
+ double a = (del2[1])*(del3[2]) - (del2[2])*(del3[1]);
+ double b = (del2[2])*(del3[0]) - (del2[0])*(del3[2]);
+ double c = (del2[0])*(del3[1]) - (del2[1])*(del3[0]);
+ int signum = sgn(a*(del1[0]) + b*(del1[1]) + c*(del1[2]));
+
+ if (fabs(signum)!= 1)
+ error->all(FLERR,"Wrong orientation in cluster of size 5"
+ "in fix filter/corotate!");
+ clist_q0[i][8] *= signum;
+ clist_q0[i][11] *= signum;
+ }
+ else
+ {
+ error->all(FLERR,"Fix filter/corotate cluster with size > 5"
+ "not yet configured...");
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * setup, needs to patch missing setup_post_force_respa() in respa...
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::setup(int vflag)
+{
+ ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+ post_force_respa(vflag,nlevels_respa-1,0);
+ ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+}
+
+void FixFilterCorotate::setup_pre_force_respa(int vflag,int ilevel){
+ pre_force_respa(vflag,ilevel,0);
+}
+
+//void FixFilterCorotate::setup_post_force_respa(int vflag,int ilevel){
+// post_force_respa(vflag,ilevel,0);
+//}
+
+double FixFilterCorotate::compute_array(int,int)
+{
+ filter_inner();
+ return 1;
+}
+
+void FixFilterCorotate::pre_force_respa(int vflag, int ilevel, int iloop)
+{
+ if (ilevel == nlevels_respa-1)
+ {
+ filter_inner();
+
+ //%Switch pointers x<-->x_filter
+ x_store = atom->x;
+ atom->x = array_atom;
+ }
+}
+
+void FixFilterCorotate::post_force_respa(int vflag, int ilevel, int iloop)
+{
+ if (ilevel == nlevels_respa-1)
+ {
+ atom->x = x_store;
+
+ comm->forward_comm_fix(this); //comm forward of forces for outer filter
+
+ filter_outer();
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixFilterCorotate::filter_inner()
+{
+ int nall = atom->nlocal + atom->nghost;
+ double **x = atom->x;
+
+ //set array_atom to atom->x, in case some atoms are not part of any cluster:
+ for (int i=0; i<nall; i++)
+ {
+ array_atom[i][0] = x[i][0];
+ array_atom[i][1] = x[i][1];
+ array_atom[i][2] = x[i][2];
+ }
+
+ int m;
+ //iterate through list:
+ for (int i=0; i<nlist; i++)
+ {
+ m = list[i];
+
+ //filter:
+ general_cluster(m,i);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixFilterCorotate::filter_outer()
+{
+ double sum1,sum2,sum3;
+ double** f = atom->f;
+ int m;
+
+ for (int i=0; i<nlist; i++)
+ {
+ m = list[i];
+ int dim = shake_flag[m];
+ if (dim == 1) //case: angle
+ dim = 3;
+
+ for (int j = 0; j < dim; j++)
+ {
+ sum1 = sum2 = sum3 = 0;
+ for (int k = 0; k < dim; k++)
+ {
+ //get local tag of k-th atom in cluster i:
+ int kk = atom->map(shake_atom[m][k]);
+ //apply derivative times force
+ sum1 += clist_derv[i][3*j][3*k]*f[kk][0]+clist_derv[i][3*j][3*k+1]*
+ f[kk][1]+clist_derv[i][3*j][3*k+2]*f[kk][2];
+ sum2 += clist_derv[i][3*j+1][3*k]*f[kk][0]+clist_derv[i][3*j+1][3*k+1]*
+ f[kk][1]+clist_derv[i][3*j+1][3*k+2]*f[kk][2];
+ sum3 += clist_derv[i][3*j+2][3*k]*f[kk][0]+clist_derv[i][3*j+2][3*k+1]*
+ f[kk][1]+clist_derv[i][3*j+2][3*k+2]*f[kk][2];
+ }
+
+ g[3*j] = sum1;
+ g[3*j+1] = sum2;
+ g[3*j+2] = sum3;
+ }
+
+ //replace force f with filtered value:
+ for (int j = 0; j < dim; j++)
+ {
+ int kk = atom->map(shake_atom[m][j]);
+
+ f[kk][0] = g[3*j];
+ f[kk][1] = g[3*j+1];
+ f[kk][2] = g[3*j+2];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * 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 FixFilterCorotate::find_clusters()
+{
+ int i,j,m,n;
+ int flag,flag_all,nbuf,size;
+ double massone;
+ tagint *buf;
+
+ if (me == 0 && screen)
+ fprintf(screen,"Finding filter clusters ...\n");
+
+ tagint *tag = atom->tag;
+ int *type = atom->type;
+ int *mask = atom->mask;
+ double *mass = atom->mass;
+ double *rmass = atom->rmass;
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+
+ int nlocal = atom->nlocal;
+ int angles_allow = atom->avec->angles_allow;
+
+ // 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_massflag[i][] = 1 if partner meets mass criterion, 0 if not
+ // 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;
+ memory->create(npartner,nlocal,"shake:npartner");
+ memory->create(nshake,nlocal,"shake:nshake");
+
+ tagint **partner_tag;
+ int **partner_mask,**partner_type,**partner_massflag;
+ int **partner_bondtype,**partner_shake,**partner_nshake;
+ memory->create(partner_tag,nlocal,max,"shake:partner_tag");
+ memory->create(partner_mask,nlocal,max,"shake:partner_mask");
+ memory->create(partner_type,nlocal,max,"shake:partner_type");
+ memory->create(partner_massflag,nlocal,max,"shake:partner_massflag");
+ memory->create(partner_bondtype,nlocal,max,"shake:partner_bondtype");
+ memory->create(partner_shake,nlocal,max,"shake:partner_shake");
+ memory->create(partner_nshake,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_massflag, partner_bondtype
+ // for bonded partners
+ // requires communication for off-proc partners
+ // -----------------------------------------------------
+
+ // fill in mask, type, massflag, bondtype if own bond partner
+ // info to store in buf for each off-proc bond = nper = 6
+ // 2 atoms IDs in bond, space for mask, type, massflag, bondtype
+ // nbufmax = largest buffer needed to hold info from any proc
+
+ int nper = 6;
+
+ 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_massflag[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];
+ if (nmass) {
+ if (rmass) massone = rmass[m];
+ else massone = mass[type[m]];
+ partner_massflag[i][j] = masscheck(massone);
+ }
+ n = bondtype_findset(i,tag[i],partner_tag[i][j],0);
+ if (n) partner_bondtype[i][j] = n;
+ else {
+ n = bondtype_findset(m,tag[i],partner_tag[i][j],0);
+ if (n) partner_bondtype[i][j] = n;
+ }
+ } else nbuf += nper;
+ }
+ }
+
+ memory->create(buf,nbuf,"filter/corotate:buf");
+
+ // 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;
+ buf[size+4] = 0;
+ n = bondtype_findset(i,tag[i],partner_tag[i][j],0);
+ if (n) buf[size+5] = n;
+ else buf[size+5] = 0;
+ size += nper;
+ }
+ }
+ }
+
+ // cycle buffer around ring of procs back to self
+
+ fsptr = this;
+ comm->ring(size,sizeof(tagint),buf,1,ring_bonds,buf);
+
+ // 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_massflag[i][j] = buf[m+4];
+ partner_bondtype[i][j] = buf[m+5];
+ m += nper;
+ }
+
+ memory->destroy(buf);
+
+ // 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(FLERR,"Did not find fix filter/corotate 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
+ // -----------------------------------------------------
+
+ int np;
+
+ for (i = 0; i < nlocal; i++) {
+ nshake[i] = 0;
+ np = npartner[i];
+ for (j = 0; j < np; 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) {
+ if (partner_massflag[i][j]) {
+ partner_shake[i][j] = 1;
+ nshake[i]++;
+ continue;
+ } else {
+ if (rmass) massone = rmass[i];
+ else massone = mass[type[i]];
+ if (masscheck(massone)) {
+ partner_shake[i][j] = 1;
+ nshake[i]++;
+ continue;
+ }
+ }
+ }
+ }
+ }
+
+ // -----------------------------------------------------
+ // 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;
+ }
+ }
+
+ memory->create(buf,nbuf,"filter/corotate:buf");
+
+ // 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
+
+ fsptr = this;
+ comm->ring(size,sizeof(tagint),buf,2,ring_nshake,buf);
+
+ // 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;
+ }
+
+ memory->destroy(buf);
+
+ // -----------------------------------------------------
+ // 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] > 4) flag = 1;
+ MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
+ if (flag_all) error->all(FLERR,"Cluster of more than 5 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(FLERR,"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_atom[i][4] = 0;
+ shake_type[i][0] = 0;
+ shake_type[i][1] = 0;
+ shake_type[i][2] = 0;
+ shake_type[i][3] = 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 && angles_allow) {
+ n = angletype_findset(i,shake_atom[i][1],shake_atom[i][2],0);
+ if (n <= 0) continue;
+ //if (angle_flag[n]) { Take all angles!
+ shake_flag[i] = 1;
+ shake_type[i][2] = 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_atom[m][4] = shake_atom[i][4];
+ shake_type[m][0] = shake_type[i][0];
+ shake_type[m][1] = shake_type[i][1];
+ shake_type[m][2] = shake_type[i][2];
+ shake_type[m][3] = shake_type[i][3];
+ } else nbuf += 11;
+ }
+ }
+
+ memory->create(buf,nbuf,"filter/corotate:buf");
+
+ // 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_atom[i][4];
+ buf[size+7] = shake_type[i][0];
+ buf[size+8] = shake_type[i][1];
+ buf[size+9] = shake_type[i][2];
+ buf[size+10] = shake_type[i][3];
+ size += 11;
+ }
+ }
+ }
+
+ // cycle buffer around ring of procs back to self
+
+ fsptr = this;
+ comm->ring(size,sizeof(tagint),buf,3,ring_shake,NULL);
+
+ memory->destroy(buf);
+
+ // -----------------------------------------------------
+ // free local memory
+ // -----------------------------------------------------
+
+ memory->destroy(npartner);
+ memory->destroy(nshake);
+ memory->destroy(partner_tag);
+ memory->destroy(partner_mask);
+ memory->destroy(partner_type);
+ memory->destroy(partner_massflag);
+ memory->destroy(partner_bondtype);
+ memory->destroy(partner_shake);
+ memory->destroy(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) {
+ * bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
+ * bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
+ * angletype_findset(i,shake_atom[i][1],shake_atom[i][2],-1);
+} else if (shake_flag[i] == 2) {
+ bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
+} else if (shake_flag[i] == 3) {
+ bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
+ bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
+} else if (shake_flag[i] == 4) {
+ bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
+ bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
+ bondtype_findset(i,shake_atom[i][0],shake_atom[i][3],-1);
+}
+}*/
+
+ // -----------------------------------------------------
+ // print info on SHAKE clusters
+ // -----------------------------------------------------
+
+ int count1,count2,count3,count4,count5;
+ count1 = count2 = count3 = count4 = count5 = 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++;
+ else if (shake_flag[i] == 5) count5++;
+ }
+
+ 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);
+ tmp = count5;
+ MPI_Allreduce(&tmp,&count5,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 size 5 clusters\n",count5/5);
+ 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 size 5 clusters\n",count5/5);
+ fprintf(logfile," %d = # of frozen angles\n",count1/3);
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * when receive buffer, scan bond partner IDs for atoms I own
+ * if I own partner:
+ * fill in mask and type and massflag
+ * search for bond with 1st atom and fill in bondtype
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::ring_bonds(int ndatum, char *cbuf)
+{
+ Atom *atom = fsptr->atom;
+ double *rmass = atom->rmass;
+ double *mass = atom->mass;
+ int *mask = atom->mask;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int nmass = fsptr->nmass;
+
+ tagint *buf = (tagint *) cbuf;
+ int m,n;
+ double massone;
+
+ for (int i = 0; i < ndatum; i += 6) {
+ m = atom->map(buf[i+1]);
+ if (m >= 0 && m < nlocal) {
+ buf[i+2] = mask[m];
+ buf[i+3] = type[m];
+ if (nmass) {
+ if (rmass) massone = rmass[m];
+ else massone = mass[type[m]];
+ buf[i+4] = fsptr->masscheck(massone);
+ }
+ if (buf[i+5] == 0) {
+ n = fsptr->bondtype_findset(m,buf[i],buf[i+1],0);
+ if (n) buf[i+5] = n;
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * when receive buffer, scan bond partner IDs for atoms I own
+ * if I own partner, fill in nshake value
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::ring_nshake(int ndatum, char *cbuf)
+{
+ Atom *atom = fsptr->atom;
+ int nlocal = atom->nlocal;
+
+ int *nshake = fsptr->nshake;
+
+ tagint *buf = (tagint *) cbuf;
+ int m;
+
+ for (int i = 0; i < ndatum; i += 3) {
+ m = atom->map(buf[i+1]);
+ if (m >= 0 && m < nlocal) buf[i+2] = nshake[m];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * when receive buffer, scan bond partner IDs for atoms I own
+ * if I own partner, fill in nshake value
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::ring_shake(int ndatum, char *cbuf)
+{
+ Atom *atom = fsptr->atom;
+ int nlocal = atom->nlocal;
+
+ int *shake_flag = fsptr->shake_flag;
+ tagint **shake_atom = fsptr->shake_atom;
+ int **shake_type = fsptr->shake_type;
+
+ tagint *buf = (tagint *) cbuf;
+ int m;
+
+ for (int i = 0; i < ndatum; i += 11) {
+ 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_atom[m][4] = buf[i+6];
+ shake_type[m][0] = buf[i+7];
+ shake_type[m][1] = buf[i+8];
+ shake_type[m][2] = buf[i+9];
+ shake_type[m][3] = buf[i+10];
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * check if massone is within MASSDELTA of any mass in mass_list
+ * return 1 if yes, 0 if not
+ * ------------------------------------------------------------------------- */
+
+int FixFilterCorotate::masscheck(double massone)
+{
+ for (int i = 0; i < nmass; i++){
+ if (fabs(mass_list[i]-massone) <= MASSDELTA) return 1;
+ }
+ return 0;
+}
+
+/* --------------------------------------------------------------------------
+ * filter one cluster of arbitrary size
+ * ---------------------------------------------------------------------------*/
+
+void FixFilterCorotate::general_cluster(int index, int index_in_list)
+{
+ //index: shake index, index_in_list: corresponding index in list
+ //get q0, nselect:
+ double *q0 = clist_q0[index_in_list];
+ int nselect1 = clist_nselect1[index_in_list];
+ int nselect2 = clist_nselect2[index_in_list];
+ int *select1 = clist_select1[index_in_list];
+ int *select2 = clist_select2[index_in_list];
+
+ int N = shake_flag[index]; //number of atoms in cluster
+ if (N == 1) N = 3; //angle cluster
+
+ double**x = atom->x;
+ double norm1, norm2, norm3;
+
+ int* list_cluster = new int[N]; // contains local IDs of cluster atoms,
+ // 0 = center
+ double* m = new double[N]; //contains local mass
+ double *r = new double[N]; //contains r[i] = 1/||del[i]||
+ double** del = new double*[N]; //contains del[i] = x_i-x_0
+ for (int i = 0; i<N; i++)
+ del[i] = new double[3];
+
+ for (int i = 0; i < N; i++)
+ {
+ list_cluster[i] = atom->map(shake_atom[index][i]);
+ m[i] = atom->mass[atom->type[list_cluster[i]]];
+ }
+
+ //%CALC r_i:
+ for (int i = 1; i < N; i++)
+ {
+ del[i][0] = x[list_cluster[i]][0] - x[list_cluster[0]][0];
+ del[i][1] = x[list_cluster[i]][1] - x[list_cluster[0]][1];
+ del[i][2] = x[list_cluster[i]][2] - x[list_cluster[0]][2];
+ domain->minimum_image(del[i]);
+ r[i] = 1.0/sqrt(del[i][0]*del[i][0]+del[i][1]*del[i][1]+
+ del[i][2]*del[i][2]);
+ }
+ //special case i=0: need del[0] for compact notation of array_atom later...
+ del[0][0] = del[0][1] = del[0][2] = r[0] = 0.0;
+
+ //calc n1:
+ n1[0] = n1[1] = n1[2] = 0.0;
+
+ int k;
+ for (int i = 0; i < nselect1; i++)
+ {
+ k = select1[i];
+ n1[0] += del[k][0]*r[k];
+ n1[1] += del[k][1]*r[k];
+ n1[2] += del[k][2]*r[k];
+ }
+ norm1 = 1.0/sqrt(n1[0]*n1[0]+n1[1]*n1[1]+n1[2]*n1[2]);
+
+ n1[0] *= norm1;
+ n1[1] *= norm1;
+ n1[2] *= norm1;
+
+ //calc n2:
+ n2[0] = n2[1] = n2[2] = 0.0;
+
+ for (int i = 0; i < nselect2; i++)
+ {
+ k = select2[i];
+ n2[0] += del[k][0]*r[k];
+ n2[1] += del[k][1]*r[k];
+ n2[2] += del[k][2]*r[k];
+ }
+ if (!nselect2) //cluster of size 2, has only one direction
+ n2[0] = n2[1] = n2[2] = 1.0;
+
+ double alpha = n2[0]*n1[0] + n2[1]*n1[1] + n2[2]*n1[2];
+
+ n2[0] -= alpha*n1[0];
+ n2[1] -= alpha*n1[1];
+ n2[2] -= alpha*n1[2];
+
+ norm2 = 1.0/sqrt(n2[0]*n2[0]+n2[1]*n2[1]+n2[2]*n2[2]);
+
+ n2[0] *= norm2;
+ n2[1] *= norm2;
+ n2[2] *= norm2;
+
+ n3[0] = n1[1]*n2[2] - n1[2]*n2[1];
+ n3[1] = n1[2]*n2[0] - n1[0]*n2[2];
+ n3[2] = n1[0]*n2[1] - n1[1]*n2[0];
+
+ norm3 = 1.0/sqrt(n3[0]*n3[0]+n3[1]*n3[1]+n3[2]*n3[2]);
+
+ n3[0] *= norm3;
+ n3[1] *= norm3;
+ n3[2] *= norm3;
+
+ //%x_filter:
+ for (int i = 0; i < N; i++)
+ {
+ k = list_cluster[i];
+ array_atom[k][0] = x[k][0]+q0[3*i]*n1[0]+q0[3*i+1]*n2[0]+q0[3*i+2]*n3[0];
+ array_atom[k][1] = x[k][1]+q0[3*i]*n1[1]+q0[3*i+1]*n2[1]+q0[3*i+2]*n3[1];
+ array_atom[k][2] = x[k][2]+q0[3*i]*n1[2]+q0[3*i+1]*n2[2]+q0[3*i+2]*n3[2];
+ }
+
+ double m_all = 0.0;
+ for (int i = 0; i<N; i++)
+ m_all += m[i];
+
+ //%x+X_center
+
+ for (int i = 0; i<N; i++)
+ {
+ k = list_cluster[i];
+ for (int j = 0; j < N; j++)
+ {
+ array_atom[k][0] += m[j]/m_all*(del[j][0]-del[i][0]);
+ array_atom[k][1] += m[j]/m_all*(del[j][1]-del[i][1]);
+ array_atom[k][2] += m[j]/m_all*(del[j][2]-del[i][2]);
+ }
+ }
+
+ //derivative:
+ //dn1dx:
+ double sum1[3][3*N];
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3*N; j++)
+ sum1[i][j] = 0;
+
+ double I3mn1n1T[3][3]; //(I_3 - n1n1T)
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3; j++)
+ I3mn1n1T[i][j] = -n1[i]*n1[j];
+ I3mn1n1T[i][i] += 1.0;
+ }
+ // sum1 part of dn1dx:
+ for (int l = 0; l < nselect1; l++)
+ {
+ k = select1[l];
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3; j++)
+ {
+ double help = del[k][i]*del[k][j]*r[k]*r[k]*r[k];
+ sum1[i][j+3*k] -= help;
+ sum1[i][j] += help;
+ }
+ sum1[i][i+3*k] += r[k];
+ sum1[i][i] -= r[k];
+ }
+ }
+ //dn1dx = norm1 * I3mn1n1T * sum1
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3*N; j++)
+ {
+ double sum = 0;
+ for (int l = 0; l<3; l++)
+ sum += I3mn1n1T[i][l]*sum1[l][j];
+ dn1dx[i][j] = norm1*sum;
+ }
+ }
+
+ //dn2dx: norm2 * I3mn2n2T * (I3mn1n1T*sum2 - rkn1pn1rk*dn1dx)
+ double sum2[3][3*N];
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3*N; j++)
+ sum2[i][j] = 0;
+
+ double I3mn2n2T[3][3]; //(I_3 - n2n2T)
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3; j++)
+ I3mn2n2T[i][j] = -n2[i]*n2[j];
+ I3mn2n2T[i][i] += 1.0;
+ }
+
+ // sum2 part of dn1dx:
+ for (int l = 0; l < nselect2; l++)
+ {
+ k = select2[l];
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3; j++)
+ {
+ double help = del[k][i]*del[k][j]*r[k]*r[k]*r[k];
+ sum2[i][j+3*k] -= help;
+ sum2[i][j] += help;
+ }
+ sum2[i][i+3*k] += r[k];
+ sum2[i][i] -= r[k];
+ }
+ }
+
+ //prefaktor:
+ double rkn1pn1rk[3][3];
+ double rk[3]; rk[0] = rk[1] = rk[2] = 0.0;
+
+ for (int i = 0; i < nselect2; i++)
+ {
+ k = select2[i];
+ rk[0] += del[k][0]*r[k];
+ rk[1] += del[k][1]*r[k];
+ rk[2] += del[k][2]*r[k];
+ }
+ //rkn1pn1rk = rkT*n1*I3 + n1*rkT
+ double scalar = rk[0]*n1[0]+rk[1]*n1[1]+rk[2]*n1[2];
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3; j++)
+ rkn1pn1rk[i][j] = n1[i]*rk[j];
+ rkn1pn1rk[i][i] += scalar;
+ }
+ //dn2dx: norm2 * I3mn2n2T * (I3mn1n1T*sum2 - rkn1pn1rk*dn1dx)
+ //sum3 = (I3mn1n1T*sum2 - rkn1pn1rk*dn1dx)
+ double sum3[3][3*N];
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3*N; j++)
+ {
+ double sum = 0;
+ for (int l = 0; l<3; l++)
+ sum += I3mn1n1T[i][l]*sum2[l][j] - rkn1pn1rk[i][l]*dn1dx[l][j];
+ sum3[i][j] = sum;
+ }
+ }
+ //dn2dx = norm2 * I3mn2n2T * sum3
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3*N; j++)
+ {
+ double sum = 0;
+ for (int l = 0; l<3; l++)
+ sum += I3mn2n2T[i][l]*sum3[l][j];
+ dn2dx[i][j] = norm2*sum;
+ }
+ }
+
+ //dn3dx = norm3 * I3mn3n3T * cross
+ double I3mn3n3T[3][3]; //(I_3 - n3n3T)
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3; j++)
+ I3mn3n3T[i][j] = -n3[i]*n3[j];
+ I3mn3n3T[i][i] += 1.0;
+ }
+
+ double cross[3][3*N];
+
+ for (int j=0; j<3*N; j++)
+ {
+ cross[0][j] = dn1dx[1][j]*n2[2] -dn1dx[2][j]*n2[1] +
+ n1[1]*dn2dx[2][j]-n1[2]*dn2dx[1][j];
+ cross[1][j] = dn1dx[2][j]*n2[0] -dn1dx[0][j]*n2[2] +
+ n1[2]*dn2dx[0][j]-n1[0]*dn2dx[2][j];
+ cross[2][j] = dn1dx[0][j]*n2[1] -dn1dx[1][j]*n2[0] +
+ n1[0]*dn2dx[1][j]-n1[1]*dn2dx[0][j];
+ }
+
+ for (int i=0; i<3; i++)
+ {
+ for (int j=0; j<3*N; j++)
+ {
+ double sum = 0;
+ for (int l = 0; l<3; l++)
+ sum += I3mn3n3T[i][l]*cross[l][j];
+ dn3dx[i][j] = norm3*sum;
+ }
+ }
+
+ for (int l=0; l<N; l++)
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3*N; j++)
+ help2[i+3*l][j] = q0[3*l]*dn1dx[i][j] + q0[3*l+1]*dn2dx[i][j] +
+ q0[3*l+2]*dn3dx[i][j];
+
+ for (int j=0; j<N; j++)
+ for (int l=0; l<N; l++)
+ for (int i=0; i<3; i++)
+ help2[i+3*l][i+3*j] += m[j]/m_all;
+
+ //need transposed of help2:
+ for (int ii = 0; ii < 3*N; ii++)
+ for (int jj = 0; jj < 3*N; jj++)
+ clist_derv[index_in_list][ii][jj] = help2[jj][ii];
+
+ //free memory
+ delete [] list_cluster;
+ delete [] m;
+ delete [] r;
+ for (int i = 0; i<N; i++)
+ delete [] del[i];
+ delete [] del;
+}
+
+int FixFilterCorotate::pack_forward_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,m;
+ double**f = atom->f;
+ m = 0;
+ for (i = 0; i < n; i++)
+ {
+ j = list[i];
+
+ buf[m++] = f[j][0];
+ buf[m++] = f[j][1];
+ buf[m++] = f[j][2];
+ }
+ return m;
+}
+
+void FixFilterCorotate::unpack_forward_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+ double **f = atom->f;
+ m = 0;
+ last = first + n;
+
+ for (i = first; i < last; i++)
+ {
+ f[i][0] = buf[m++];
+ f[i][1] = buf[m++];
+ f[i][2] = buf[m++];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * find a bond between global atom IDs n1 and n2 stored with local atom i
+ * if find it:
+ * if setflag = 0, return bond type
+ * if setflag = -1/1, set bond type to negative/positive and return 0
+ * if do not find it, return 0
+ * ------------------------------------------------------------------------- */
+
+int FixFilterCorotate::bondtype_findset(int i, tagint n1, tagint n2,
+ int setflag)
+{
+ int m,nbonds;
+
+ tagint *tag = atom->tag;
+ tagint **bond_atom = atom->bond_atom;
+ nbonds = atom->num_bond[i];
+
+ 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) {
+ if (setflag == 0)
+ return atom->bond_type[i][m];
+
+ if ((setflag < 0 && atom->bond_type[i][m] > 0) ||
+ (setflag > 0 && atom->bond_type[i][m] < 0))
+ atom->bond_type[i][m] = -atom->bond_type[i][m];
+
+ }
+
+ return 0;
+}
+
+/* ----------------------------------------------------------------------
+ * find an angle with global end atom IDs n1 and n2 stored with local atom i
+ * if find it:
+ * if setflag = 0, return angle type
+ * if setflag = -1/1, set angle type to negative/positive and return 0
+ * if do not find it, return 0
+ * ------------------------------------------------------------------------- */
+
+int FixFilterCorotate::angletype_findset(int i, tagint n1, tagint n2,
+ int setflag)
+{
+ int m,nangles;
+
+ tagint **angle_atom1 = atom->angle_atom1;
+ tagint **angle_atom3 = atom->angle_atom3;
+ nangles = atom->num_angle[i];
+
+ 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) {
+ if (setflag == 0) {
+ return atom->angle_type[i][m];
+ }
+
+ if ((setflag < 0 && atom->angle_type[i][m] > 0) ||
+ (setflag > 0 && atom->angle_type[i][m] < 0))
+ atom->angle_type[i][m] = -atom->angle_type[i][m];
+ }
+
+ return 0;
+}
+
+/* ----------------------------------------------------------------------
+ * allocate local atom-based arrays
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::grow_arrays(int nmax)
+{
+ memory->grow(array_atom,nmax,3,"FilterCorotate:peratomarray");
+ memory->grow(shake_flag,nmax,"FilterCorotate::shake_flag");
+ memory->grow(shake_atom,nmax,5,"FilterCorotate::shake_atom");
+ memory->grow(shake_type,nmax,4,"FilterCorotate::shake_type");
+}
+
+/* ----------------------------------------------------------------------
+ * memory usage of local atom-based arrays
+ * ------------------------------------------------------------------------- */
+
+double FixFilterCorotate::memory_usage()
+{
+ double bytes = 0;
+ //GROW:
+ bytes += 3*sizeof(double) + 5*sizeof(tagint) + 5*sizeof(int);
+ //clist
+ bytes += 13*atom->nlocal*sizeof(int);
+ bytes += 15*16*nlocal*sizeof(double);
+
+ //fixed:
+ int nb = atom->nbondtypes+1;
+ int na = atom->nangletypes+1;
+ int nt = atom->ntypes+1;
+ bytes += (nb+na+nt)*sizeof(int);
+ bytes += (nt-1+nb+na+15*15+18+10*15)*sizeof(double);
+
+ //output:
+ //bytes += (2*nb+2*na)*sizeof(int);
+ //bytes += (6*na+6*nb)*sizeof(double);
+
+ return bytes;
+}
+
+
+/* ----------------------------------------------------------------------
+ * copy values within local atom-based arrays
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::copy_arrays(int i, int j, int delflag)
+{
+ 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];
+ shake_type[j][2] = shake_type[i][2];
+ } 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];
+ } else if (flag == 5) {
+ 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_atom[j][4] = shake_atom[i][4];
+ shake_type[j][0] = shake_type[i][0];
+ shake_type[j][1] = shake_type[i][1];
+ shake_type[j][2] = shake_type[i][2];
+ shake_type[j][3] = shake_type[i][3];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ * initialize one atom's array values, called when atom is created
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::set_arrays(int i)
+{
+ shake_flag[i] = 0;
+}
+
+/* ----------------------------------------------------------------------
+ * update one atom's array values
+ * called when molecule is created from fix gcmc
+ * ------------------------------------------------------------------------- */
+
+void FixFilterCorotate::update_arrays(int i, int atom_offset)
+{
+ int flag = shake_flag[i];
+
+ if (flag == 1) {
+ shake_atom[i][0] += atom_offset;
+ shake_atom[i][1] += atom_offset;
+ shake_atom[i][2] += atom_offset;
+ } else if (flag == 2) {
+ shake_atom[i][0] += atom_offset;
+ shake_atom[i][1] += atom_offset;
+ } else if (flag == 3) {
+ shake_atom[i][0] += atom_offset;
+ shake_atom[i][1] += atom_offset;
+ shake_atom[i][2] += atom_offset;
+ } else if (flag == 4) {
+ shake_atom[i][0] += atom_offset;
+ shake_atom[i][1] += atom_offset;
+ shake_atom[i][2] += atom_offset;
+ shake_atom[i][3] += atom_offset;
+ } else if (flag == 5) {
+ shake_atom[i][0] += atom_offset;
+ shake_atom[i][1] += atom_offset;
+ shake_atom[i][2] += atom_offset;
+ shake_atom[i][3] += atom_offset;
+ shake_atom[i][4] += atom_offset;
+ }
+}
+
+
+/* ----------------------------------------------------------------------
+ * pack values in local atom-based arrays for exchange with another proc
+ * ------------------------------------------------------------------------- */
+
+int FixFilterCorotate::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];
+ buf[m++] = shake_type[i][2];
+ } 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];
+ } else if (flag == 5) {
+ 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_atom[i][4];
+ buf[m++] = shake_type[i][0];
+ buf[m++] = shake_type[i][1];
+ buf[m++] = shake_type[i][2];
+ buf[m++] = shake_type[i][3];
+ }
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ * unpack values in local atom-based arrays from exchange with another proc
+ * ------------------------------------------------------------------------- */
+
+int FixFilterCorotate::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<tagint> (buf[m++]);
+ shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][2] = static_cast<tagint> (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<tagint> (buf[m++]);
+ shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
+ shake_type[nlocal][0] = static_cast<int> (buf[m++]);
+ } else if (flag == 3) {
+ shake_atom[nlocal][0] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][2] = static_cast<tagint> (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 == 4) {
+ shake_atom[nlocal][0] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][2] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][3] = static_cast<tagint> (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 == 5) {
+ shake_atom[nlocal][0] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][1] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][2] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][3] = static_cast<tagint> (buf[m++]);
+ shake_atom[nlocal][4] = static_cast<tagint> (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++]);
+ shake_type[nlocal][3] = static_cast<int> (buf[m++]);
+ }
+ return m;
+}
\ No newline at end of file
diff --git a/src/USER-MISC/fix_filter_corotate.h b/src/USER-MISC/fix_filter_corotate.h
index fdac4bed748374bfaf97ef97f3476211cde5b8b0..da9e2dedb3eb19c421201e5c00c8c833b33a1df9 100644
--- a/src/USER-MISC/fix_filter_corotate.h
+++ b/src/USER-MISC/fix_filter_corotate.h
@@ -1,139 +1,139 @@
-/* ----------------------------------------------------------------------
- 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.
- ------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
- Contributing author: Lukas Fath (KIT)
- some subroutines are from fix_shake.cpp
- ------------------------------------------------------------------------- */
-
-#ifdef FIX_CLASS
-
-FixStyle(filter/corotate,FixFilterCorotate)
-
-#else
-
-#ifndef LMP_FIX_FILTER_COROTATE_H
-#define LMP_FIX_FILTER_COROTATE_H
-
-#include "fix.h"
-
-namespace LAMMPS_NS
-{
-
- class FixFilterCorotate : public Fix
- {
- public:
-
- FixFilterCorotate(class LAMMPS *, int, char **);
- ~FixFilterCorotate();
- void setup(int);
- void setup_pre_neighbor();
- void pre_neighbor();
- void setup_pre_force_respa(int,int);
-// void setup_post_force_respa(int,int);
- void pre_force_respa(int, int, int);
- void post_force_respa(int, int, int);
-
- void init();
- int setmask();
-
- double compute_array(int,int);
-
- int pack_forward_comm(int, int *, double *, int, int *);
- void unpack_forward_comm(int, int, double *);
- void grow_arrays(int );
- double memory_usage();
-
- void copy_arrays(int, int, int);
- void set_arrays(int);
- void update_arrays(int, int);
-
- int pack_exchange(int, double *);
- int unpack_exchange(int, double *);
-
- protected:
-
- int me,nprocs;
-
- int flevel; //filtered respa level
-
- double **help2; //temp derivative
- double **x_store; //temp for atom->x
- double *g; //temp for derivative
-
- double*n1,*n2,*n3, *del1, *del2,*del3;
-
- double**dn1dx,**dn2dx,**dn3dx;
-
- int *bond_flag,*angle_flag; // bond/angle types to constrain
- int *type_flag; // constrain bonds to these types
- double *mass_list; // constrain bonds to these masses
- int nmass; // # of masses in mass_list
-
- int molecular; // copy of atom->molecular
- double *bond_distance,*angle_distance; // constraint distances
-
- int nlevels_respa; // copies of needed rRESPA variables
-
- double **x,**v,**f; // local ptrs to atom class quantities
- double *mass,*rmass;
- int *type;
- int nlocal;
- // atom-based arrays
- int *shake_flag; // 0 if atom not in SHAKE cluster
- // 1 = size 3 angle cluster
- // 2,3,4 = size of bond-only cluster
- tagint **shake_atom; // global IDs of atoms in cluster
- // central atom is 1st
- // lowest global ID is 1st for size 2
- int **shake_type; // bondtype of each bond in cluster
- // for angle cluster, 3rd value
- // is angletype
- int *nshake; // count
-
- int *list; // list of clusters to SHAKE
- int nlist,maxlist; // size and max-size of list
- double ***clist_derv; //stores derivative
- double **clist_q0; //stores reference config
- int *clist_nselect1, *clist_nselect2; //stores length of each selec. list
- int **clist_select1, **clist_select2; //stores selection lists
-
- void find_clusters();
- int masscheck(double);
-
- void filter_inner();
- void filter_outer();
-
- void general_cluster(int,int);
-
- void stats();
- int bondtype_findset(int, tagint, tagint, int);
- int angletype_findset(int, tagint, tagint, int);
-
- // static variable for ring communication callback to access class data
- // callback functions for ring communication
-
- static FixFilterCorotate *fsptr;
- static void ring_bonds(int, char *);
- static void ring_nshake(int, char *);
- static void ring_shake(int, char *);
-
- int sgn(double val) {
- return (0 < val) - (val < 0);
- };
- };
-
-}
-
-#endif
+/* ----------------------------------------------------------------------
+ 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.
+ ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Lukas Fath (KIT)
+ some subroutines are from fix_shake.cpp
+ ------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(filter/corotate,FixFilterCorotate)
+
+#else
+
+#ifndef LMP_FIX_FILTER_COROTATE_H
+#define LMP_FIX_FILTER_COROTATE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS
+{
+
+ class FixFilterCorotate : public Fix
+ {
+ public:
+
+ FixFilterCorotate(class LAMMPS *, int, char **);
+ ~FixFilterCorotate();
+ void setup(int);
+ void setup_pre_neighbor();
+ void pre_neighbor();
+ void setup_pre_force_respa(int,int);
+// void setup_post_force_respa(int,int);
+ void pre_force_respa(int, int, int);
+ void post_force_respa(int, int, int);
+
+ void init();
+ int setmask();
+
+ double compute_array(int,int);
+
+ int pack_forward_comm(int, int *, double *, int, int *);
+ void unpack_forward_comm(int, int, double *);
+ void grow_arrays(int );
+ double memory_usage();
+
+ void copy_arrays(int, int, int);
+ void set_arrays(int);
+ void update_arrays(int, int);
+
+ int pack_exchange(int, double *);
+ int unpack_exchange(int, double *);
+
+ protected:
+
+ int me,nprocs;
+
+ int flevel; //filtered respa level
+
+ double **help2; //temp derivative
+ double **x_store; //temp for atom->x
+ double *g; //temp for derivative
+
+ double*n1,*n2,*n3, *del1, *del2,*del3;
+
+ double**dn1dx,**dn2dx,**dn3dx;
+
+ int *bond_flag,*angle_flag; // bond/angle types to constrain
+ int *type_flag; // constrain bonds to these types
+ double *mass_list; // constrain bonds to these masses
+ int nmass; // # of masses in mass_list
+
+ int molecular; // copy of atom->molecular
+ double *bond_distance,*angle_distance; // constraint distances
+
+ int nlevels_respa; // copies of needed rRESPA variables
+
+ double **x,**v,**f; // local ptrs to atom class quantities
+ double *mass,*rmass;
+ int *type;
+ int nlocal;
+ // atom-based arrays
+ int *shake_flag; // 0 if atom not in SHAKE cluster
+ // 1 = size 3 angle cluster
+ // 2,3,4 = size of bond-only cluster
+ tagint **shake_atom; // global IDs of atoms in cluster
+ // central atom is 1st
+ // lowest global ID is 1st for size 2
+ int **shake_type; // bondtype of each bond in cluster
+ // for angle cluster, 3rd value
+ // is angletype
+ int *nshake; // count
+
+ int *list; // list of clusters to SHAKE
+ int nlist,maxlist; // size and max-size of list
+ double ***clist_derv; //stores derivative
+ double **clist_q0; //stores reference config
+ int *clist_nselect1, *clist_nselect2; //stores length of each selec. list
+ int **clist_select1, **clist_select2; //stores selection lists
+
+ void find_clusters();
+ int masscheck(double);
+
+ void filter_inner();
+ void filter_outer();
+
+ void general_cluster(int,int);
+
+ void stats();
+ int bondtype_findset(int, tagint, tagint, int);
+ int angletype_findset(int, tagint, tagint, int);
+
+ // static variable for ring communication callback to access class data
+ // callback functions for ring communication
+
+ static FixFilterCorotate *fsptr;
+ static void ring_bonds(int, char *);
+ static void ring_nshake(int, char *);
+ static void ring_shake(int, char *);
+
+ int sgn(double val) {
+ return (0 < val) - (val < 0);
+ };
+ };
+
+}
+
+#endif
#endif
\ No newline at end of file