Commit db0281b4 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

Merge branch 'filter_corotate' of https://github.com/lukin17/lammps into pull-416

parents 0d8f74f0 2f5e711a

src/USER-MISC/fix_filter_corotate.cpp

+2062 −2065
Original line number Diff line number Diff line
@@ -17,8 +17,8 @@ 
   ------------------------------------------------------------------------- */



#include <mpi.h>

#include "string.h"

#include "stdlib.h"

#include <string.h>

#include <stdlib.h>

#include "fix_filter_corotate.h"

#include "atom.h"

#include "atom_vec.h"
@@ -39,7 +39,6 @@ 
#include "domain.h"

#include "integrate.h"

#include "respa.h"

#include "timer.h"

#include "neighbor.h"

#include "citeme.h"
@@ -57,7 +56,8 @@ FixFilterCorotate *FixFilterCorotate::fsptr; 
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"

" 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"
@@ -71,7 +71,8 @@ static const char cite_filter_corotate[] = 


/* ---------------------------------------------------------------------- */



FixFilterCorotate::FixFilterCorotate(LAMMPS *lmp, int narg, char **arg):Fix(lmp,narg,arg)

FixFilterCorotate::FixFilterCorotate(LAMMPS *lmp, int narg, char **arg) :

  Fix(lmp,narg,arg)

{

  if (lmp->citeme) lmp->citeme->add(cite_filter_corotate);
@@ -80,9 +81,11 @@ FixFilterCorotate::FixFilterCorotate(LAMMPS *lmp, int narg, char **arg):Fix(lmp, 


  molecular = atom->molecular;

  if (molecular == 0)

    error->all(FLERR,"Cannot use fix filter/corotate with non-molecular system");

    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");

    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
@@ -135,7 +138,8 @@ FixFilterCorotate::FixFilterCorotate(LAMMPS *lmp, int narg, char **arg):Fix(lmp, 


    } 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 (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;
@@ -164,7 +168,6 @@ FixFilterCorotate::FixFilterCorotate(LAMMPS *lmp, int narg, char **arg):Fix(lmp, 


  x_store = NULL;



  

  //STUFF

  g = NULL;

  help2 = NULL;
@@ -267,15 +270,18 @@ void FixFilterCorotate::init() 
  int i;

  // error if more than one filter

  int count = 0;

  for (i = 0; i < modify->nfix; i++)

  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++)

  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 (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
@@ -357,7 +363,7 @@ void FixFilterCorotate::pre_neighbor() 


  nlist = 0;



  for (int i = 0; i < nlocal; i++)

  for (int i = 0; i < nlocal; i++){

    if (shake_flag[i]) {

      if (shake_flag[i] == 2) {

        atom1 = atom->map(shake_atom[i][0]);
@@ -408,7 +414,8 @@ void FixFilterCorotate::pre_neighbor() 
        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) {

        if (atom1 == -1 || atom2 == -1 || atom3 == -1 ||

          atom4 == -1 || atom5 == -1) {

          char str[128];

          sprintf(str,"Cluster atoms "

          TAGINT_FORMAT " " TAGINT_FORMAT " "
@@ -423,10 +430,11 @@ void FixFilterCorotate::pre_neighbor() 
          list[nlist++] = i;

      }

    }

  }



  //also build reference config lists

  int m,N,oxy;

    double r0,r1,r2,r3,theta0;

  double r0,r1,r2,theta0;

  double m1,m2,m3,m4,m5,m_all;



  double *mass = atom->mass;
@@ -442,7 +450,8 @@ void FixFilterCorotate::pre_neighbor() 
    {

      //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);

        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
@@ -587,7 +596,9 @@ void FixFilterCorotate::pre_neighbor() 
      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!");

      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;
@@ -659,15 +670,17 @@ void FixFilterCorotate::pre_neighbor() 
      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!");

      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...");

      error->all(FLERR,"Fix filter/corotate cluster with size > 5"

        "not yet configured...");

    }

      

  }

}
@@ -682,7 +695,6 @@ void FixFilterCorotate::setup(int vflag) 
  ((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);

}
@@ -691,7 +703,6 @@ void FixFilterCorotate::setup_pre_force_respa(int vflag,int ilevel){ 
//  post_force_respa(vflag,ilevel,0);

//}





double FixFilterCorotate::compute_array(int,int)

{

  filter_inner();
@@ -710,7 +721,6 @@ void FixFilterCorotate::pre_force_respa(int vflag, int ilevel, int iloop) 
  }

}





void FixFilterCorotate::post_force_respa(int vflag, int ilevel, int iloop)

{

  if (ilevel == nlevels_respa-1)
@@ -727,15 +737,9 @@ void FixFilterCorotate::post_force_respa(int vflag, int ilevel, int iloop) 


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++)

  {
@@ -753,19 +757,14 @@ void FixFilterCorotate::filter_inner() 
    //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++)
@@ -783,9 +782,12 @@ void FixFilterCorotate::filter_outer() 
        //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];

        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;
@@ -802,12 +804,9 @@ void FixFilterCorotate::filter_outer() 
      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
@@ -818,17 +817,14 @@ void FixFilterCorotate::filter_outer() 


void FixFilterCorotate::find_clusters()

{

  int i,j,m,n,iatom;

  int i,j,m,n;

  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;
@@ -837,7 +833,6 @@ void FixFilterCorotate::find_clusters() 
  int **nspecial = atom->nspecial;

  tagint **special = atom->special;



  

  int nlocal = atom->nlocal;

  int angles_allow = atom->avec->angles_allow;
@@ -885,14 +880,12 @@ void FixFilterCorotate::find_clusters() 
  // 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
@@ -965,8 +958,9 @@ void FixFilterCorotate::find_clusters() 
  m = 0;

  while (m < size) {

    i = atom->map(buf[m]);

    for (j = 0; j < npartner[i]; j++)

    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];
@@ -985,16 +979,18 @@ void FixFilterCorotate::find_clusters() 
  // else it's an error



  flag = 0;

  for (i = 0; i < nlocal; i++)

  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");

  if (flag_all)

    error->all(FLERR,"Did not find fix filter/corotate partner info");



  // -----------------------------------------------------

  // identify SHAKEable bonds
@@ -1091,8 +1087,9 @@ void FixFilterCorotate::find_clusters() 
  m = 0;

  while (m < size) {

    i = atom->map(buf[m]);

    for (j = 0; j < npartner[i]; j++)

    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;

  }
@@ -1150,8 +1147,9 @@ void FixFilterCorotate::find_clusters() 
    shake_type[i][3] = 0;



    if (nshake[i] == 1) {

      for (j = 0; j < npartner[i]; j++)

      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];
@@ -1429,8 +1427,9 @@ void FixFilterCorotate::ring_shake(int ndatum, char *cbuf) 


int FixFilterCorotate::masscheck(double massone)

{

  for (int i = 0; i < nmass; i++)

  for (int i = 0; i < nmass; i++){

    if (fabs(mass_list[i]-massone) <= MASSDELTA) return 1;

  }

  return 0;

}
@@ -1441,7 +1440,6 @@ int FixFilterCorotate::masscheck(double massone) 
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];
@@ -1455,7 +1453,8 @@ void FixFilterCorotate::general_cluster(int index, int index_in_list) 
  double**x = atom->x;

  double norm1, norm2, norm3;



  int* list_cluster = new int[N];  //contains local IDs of cluster atoms, 0 = center

  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
@@ -1475,7 +1474,8 @@ void FixFilterCorotate::general_cluster(int index, int index_in_list) 
    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]);

    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;
@@ -1688,9 +1688,12 @@ void FixFilterCorotate::general_cluster(int index, int index_in_list) 


  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];

    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++)
@@ -1703,10 +1706,12 @@ void FixFilterCorotate::general_cluster(int index, int index_in_list) 
      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];   

        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++)
@@ -1727,9 +1732,8 @@ void FixFilterCorotate::general_cluster(int index, int index_in_list) 
  delete [] del;

}







int FixFilterCorotate::pack_forward_comm(int n, int *list, double *buf,int pbc_flag, int *pbc)

int FixFilterCorotate::pack_forward_comm(int n, int *list, double *buf,

                                         int pbc_flag, int *pbc)

{

  int i,j,m;

  double**f = atom->f;
@@ -1768,11 +1772,10 @@ void FixFilterCorotate::unpack_forward_comm(int n, int first, double *buf) 
 *  if do not find it, return 0

 * ------------------------------------------------------------------------- */



int FixFilterCorotate::bondtype_findset(int i, tagint n1, tagint n2, int setflag)

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;
@@ -1783,8 +1786,6 @@ int FixFilterCorotate::bondtype_findset(int i, tagint n1, tagint n2, int setflag 
    if (n1 == bond_atom[i][m] && n2 == tag[i]) break;

  }



  

  

  if (m < nbonds) {

    if (setflag == 0)

      return atom->bond_type[i][m];
@@ -1806,11 +1807,10 @@ int FixFilterCorotate::bondtype_findset(int i, tagint n1, tagint n2, int setflag 
 *  if do not find it, return 0

 * ------------------------------------------------------------------------- */



int FixFilterCorotate::angletype_findset(int i, tagint n1, tagint n2, int setflag)

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;
@@ -1821,7 +1821,6 @@ int FixFilterCorotate::angletype_findset(int i, tagint n1, tagint n2, int setfla 
    if (n1 == angle_atom3[i][m] && n2 == angle_atom1[i][m]) break;

  }



  

  if (m < nangles) {

    if (setflag == 0) {

      return atom->angle_type[i][m];
@@ -1830,13 +1829,11 @@ int FixFilterCorotate::angletype_findset(int i, tagint n1, tagint n2, int setfla 
    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

 * ------------------------------------------------------------------------- */

src/USER-MISC/fix_filter_corotate.h

+138 −138

File changed.

Contains only whitespace changes.

CRA Git | Maintained and supported by SUSTech CRA and CCSE