fix_topo2.cpp

/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   https://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 authors: Yair Augusto Gutiérrez Fosado
                         Davide Michieletto
------------------------------------------------------------------------- */
#include "fix_topo2.h"

//The atom class provides access to various atom properties, including the atom type, molecular ID, position, velocity, and force
#include "atom.h"

//Came by default in fix_swap_atom
#include <cmath>
#include <cctype>
#include <cfloat>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "compute.h"
#include "group.h"
#include "domain.h"
#include "region.h"
#include "random_park.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"

//Davide include them (with more)
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "random_mars.h"
#include "citeme.h"
#include <iostream>
#include<stdlib.h>
#include <time.h>       /* Important for random number generator */
#include <sstream>      // std::stringstream
#include<fstream>       //ofstream

#include "input.h"
#include "variable.h"




using namespace std;
using namespace LAMMPS_NS;
using namespace FixConst;

static const char cite_fix_topo2[] =
"Non-equilibrium Topo2 Strand Crossing:\n\n"
  "@Article{Michieletto2018,\n"
  " author = {Michieletto et al},\n"
  " title = {xx},\n"
  " journal = {xx},\n"
  " year =    2018,\n"
  " volume =  xx,\n"
  " pages =   {xx}\n"
  "}\n\n";

/* ---------------------------------------------------------------------- */
//Constructor: The name of the class is FixTopo2 and the fix-style-name to 
//be used in the lammps script is fix_topo2.
FixTopo2::FixTopo2(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg),
  alist(nullptr), list(nullptr), random_equal(nullptr), random_unequal(nullptr), random(nullptr)
{
  //Cite
  if (lmp->citeme) lmp->citeme->add(cite_fix_topo2);
  
  //Check that the atom_style is molecular
  if (atom->molecular != 1)
    error->all(FLERR,"Cannot use fix_topo2 with non-molecular systems");

  //Number of arguments for the fix. The first three arguments are parsed by Fix base class constructor.
  //The rest are specific to this fix. From 4 to 8 are mandatory for all cases.
  //4.- flagtopo: Static or randjump or jump2maxdens or jum2marcurv
  //5.- ntopo: the number of topo2 in the system.
  //6.- ltopo: the length of each topo2 in the system.
  //7.- topotype: the topoII particle type.
  //8.- seed: random seed
  
  //If Static we need one more argument
  //9.-tsteptopo: at which timestep we want to introduce the topo2
  
  //If randjump or jump2maxdens or jum2marcurv, we need three arguments
  //9.- nevery: how often the fix is called.
  //10.- prob: probability that the topo2 changes position each nevery timesteps.
  //11.- intertype: intermediate type assigned to atoms from a section where topo 2 is jumping from.
    
  //commands use in fix_bond_swap
  //force_reneighbor = -1;
  next_reneighbor = -1;
  //next_reneighbor = 1;
  vector_flag = 1;
  size_vector = 2;
  global_freq = 1;
  extvector = 0;
   

  //Required arguments: common to all cases
  flagtopo=0;
  if (strcmp(arg[3],"static") == 0) {flagtopo=1;}
  else if (strcmp(arg[3],"randjump") == 0) {flagtopo=2;}
  else if (strcmp(arg[3],"jump2maxdens") == 0) {flagtopo=3;}
  else if (strcmp(arg[3],"jump2maxcurv") == 0) {flagtopo=4;}
  if( flagtopo==0 ){error->all(FLERR,"Illegal fix topo2 command");}

  //Check that the number of arguments is the correct one
  if( flagtopo==1){
      if(narg!=9){error->all(FLERR,"Illegal fix topo2 command, a different number of arguments is required for static");}
  }
  if( flagtopo==2){
      if(narg!=11){error->all(FLERR,"Illegal fix topo2 command, a different number of arguments is required for randjump");}
  }
  if( flagtopo==3){
      if(narg!=11){error->all(FLERR,"Illegal fix topo2 command, a different number of arguments is required for jump2maxdens");}
  }
  if( flagtopo==4){
      if(narg!=11){error->all(FLERR,"Illegal fix topo2 command, a different number of arguments is required for jump2maxcurv");}
  }
  
  ntopo = utils::inumeric(FLERR,arg[4],false,lmp);
  if (ntopo < 0 || ntopo > atom->natoms) error->all(FLERR,"Illegal fix topo2 command: in ntopo");

  ltopo = utils::inumeric(FLERR,arg[5],false,lmp);
  int length = ntopo*ltopo;
  if (ltopo < 0 || length > atom->natoms) error->all(FLERR,"Illegal fix topo2 command: in ltopo");
  
  topotype = utils::inumeric(FLERR,arg[6],false,lmp);
  int ilo,ihi,jlo,jhi;
  utils::bounds(FLERR,arg[6],1,atom->ntypes,ilo,ihi,error);
  
  seed = utils::inumeric(FLERR,arg[7],false,lmp);
  if (seed <= 0) error->all(FLERR,"Illegal fix topo2 command");
  
  //In case static, we have to set the timestep at which we want to introduce the topo2
  if(flagtopo==1){
      tsteptopo = utils::inumeric(FLERR,arg[8],false,lmp);
      if (tsteptopo < 0) error->all(FLERR,"Illegal fix topo2 command");
      if (tsteptopo < update->beginstep) error->all(FLERR,"Illegal fix topo2 command");
  }
  
  //In case of the last three cases additional arguments needed:
  if(flagtopo>=2){
      nevery = utils::inumeric(FLERR,arg[8],false,lmp);
      if (nevery <= 0) error->all(FLERR,"Illegal fix topo2 command");
  
      prob = utils::numeric(FLERR,arg[9],false,lmp);
      if (prob<0.0 || prob>1.0) error->all(FLERR,"Illegal fix topo2 command");
      
      intertype = utils::inumeric(FLERR,arg[10],false,lmp);
      utils::bounds(FLERR,arg[10],1,atom->ntypes,ilo,ihi,error);
  }

   
  // random number generator, same for all procs
  random_equal = new RanPark(lmp,seed);

  // random number generator, not the same for all procs
  random_unequal = new RanPark(lmp,seed);
   
  // initialize Marsaglia RNG with processor-unique seed
  random = new RanMars(lmp,seed + comm->me);
  
  //To get a different random number every time the program is executed
  srand(time(NULL)*seed);
  
  //Decleare selection list 
  select = nullptr;
  beads  = nullptr;
  
  //Declare atom list
  alist = nullptr;
}



/****************/
/*  Destructor  */
/****************/
FixTopo2::~FixTopo2()
{
  delete random_equal;
  delete random_unequal;
  delete random;
  
  memory->destroy(select);
  memory->destroy(beads);
  memory->destroy(alist);
}



/**************/
/*  Set mask  */
/**************/
//It determines the stage in the timestep at which the fix will be executed. 
//The post integrate method is set for operations that need to happen immediately after updates by for example fix nve which performs the start-of-timestep velocity-Verlet integration operations to update velocities by a half-step, and coordinates by a full step. Only a few fixes use this, e.g. to reflect particles off box boundaries.
int FixTopo2::setmask()
{
  int mask = 0;
  mask |= POST_INTEGRATE;
  return mask;
}



/***********************/
/*  The init() method  */
/***********************/
//This is called at the beginning of each run, simply sets some internal flags
void FixTopo2::init()
{               
  // require an atom style with molecule IDs
  if (atom->molecule == NULL)
    error->all(FLERR, "Must use atom style with molecule IDs with fix topo2");

  // pair and bonds must be defined
  if (force->pair == NULL || force->bond == NULL)
    error->all(FLERR,"Fix topo2 requires pair and bond styles");
    
  // angles must be defined
  if (force->angle == NULL && atom->nangles > 0 && comm->me == 0)
  error->warning(FLERR,"Fix topo2 will ignore defined angles");

  if (force->pair->single_enable == 0)
    error->all(FLERR,"Pair style does not support fix topo2");

  // no dihedral or improper potentials allowed
  if (force->dihedral || force->improper)
    error->all(FLERR,"Fix topo2 cannot use dihedral or improper styles");

  // special bonds must be 0 1 1
  if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 ||
      force->special_lj[3] != 1.0)
    error->all(FLERR,"Fix topo2 requires special_bonds = 0,1,1");
    

/*    
  // need a full neighbor list, built every Nevery steps
  int irequest = neighbor->request(this,instance_me);
  neighbor->requests[irequest]->pair = 0;
  neighbor->requests[irequest]->fix = 1;
  neighbor->requests[irequest]->half = 0;
  neighbor->requests[irequest]->full = 1;
  //neighbor->requests[irequest]->occasional = 1;
*/
}

//Then, an init_list() method needs to be introduced to get hold of the neighbor list pointer, with a placeholder named list
//void FixTopo2::init_list(int /*id*/, NeighList *ptr)
//{
//  list = ptr;
//}

/************************************************/
/*  The action of this fix is implemented here  */
/************************************************/
void FixTopo2::post_integrate()
{   
    //Print flag (1) yes print, (0) no
    int printflag=1;
    stringstream writeFile;
    ofstream write;
    
    //Rankid
    int myrank;
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
        
    //Static: TopoII chooses a random position and it stays there until the end of the simulation
    if(flagtopo==1){
        if(tsteptopo != update->ntimestep ) return;
        
        //This is to print to a file with the selected ramdom id (for static).
        if(printflag==1 && myrank==0){
            writeFile <<"static_fromLAMMPS.dat";
            write.open(writeFile.str().c_str(), std::ios_base::app);
        }
        
        //Selects a random atom id
        int rng = randomid();
        if(printflag==1 && myrank==0){write << update->ntimestep << " " << rng << " " << topotype << endl;}
        
        //Change the bead-type of the particle with the previous id (from 1 to topotype)
        //and also to the (ltopo-1) particles that are behind.
        placetopo(rng);
    }
    
    
    //Random jump: with a certain frequency (nevery) the topoII can jump to a random position.
    //The decision of jumping is done with a certain probability.
    if(flagtopo==2){        
        if (update->ntimestep % nevery) return;
        
        //This is to print to a file with the selected ramdom id (for randjump)
        if(printflag==1 && myrank==0){
            writeFile <<"randjump_fromLAMMPS.dat";
            write.open(writeFile.str().c_str(), std::ios_base::app);
        }
        
        //The first time the programm is called it introduces a topoII with probability=1.
        if (update->ntimestep == nevery){
            int rng = randomid();
            if(printflag==1 && myrank==0){write << update->ntimestep << " " << rng << " " << topotype << endl;}
            placetopo(rng);
        }
        
        //The next times it is called
        if (update->ntimestep > nevery){
            //Probability of moving
            double pmv = pmove();
            
            if(pmv<=prob){
                int rng = randomid();
                if(printflag==1 && myrank==0){write << update->ntimestep << " " << rng << " " << topotype << endl;}
                changetype();
                placetopo(rng);               
            }
        }
    }

    
    //Jump to maximum density: with a certain frequency (nevery) the topoII can jump to the position of maximum bead density.
    //The decision of jumping is done with a certain probability.
    if(flagtopo==3){
        if (update->ntimestep % nevery) return;
        
        //This is to print to a file with the bead-id with maximum density
        if(printflag==1 && myrank==0){
            writeFile <<"jump2maxdens_fromLAMMPS.dat";
            write.open(writeFile.str().c_str(), std::ios_base::app);
        }
        
        //Window size (in beads) within which it is reasonable to accommodate a knot. Smaller values will not be informative
        int w=50;

        //The first time the programm is called it introduces a topoII with probability=1.
        if (update->ntimestep == nevery){
            int idmd = maxdensid(w);
            if(printflag==1 && myrank==0){write << update->ntimestep << " " << idmd << " " << topotype << endl;}
            placetopo(idmd);
        }

        //The next times it is called
        if (update->ntimestep > nevery){
            //Probability of moving
            double pmv = pmove();

            if(pmv<=prob){
                int idmd = maxdensid(w);
                if(printflag==1 && myrank==0){write << update->ntimestep << " " << idmd << " " << topotype << endl;}
                changetype();
                placetopo(idmd);
            }
        }
    }


    //Jump to maximum curvature: with a certain frequency (nevery) the topoII can jump to the position of maximum curvature.
    //The decision of jumping is done with a certain probability.
    if(flagtopo==4){    
        if (update->ntimestep % nevery) return;
        
        //This is to print to a file with the bead-id with maximum curvature
        if(printflag==1 && myrank==0){
            writeFile <<"jump2maxcurv_fromLAMMPS.dat";
            write.open(writeFile.str().c_str(), std::ios_base::app);
        }
        
        //Window size (in beads) within which the local curvature will be computed and integrated
        int w=50;
                
        //The first time the programm is called it introduces a topoII with probability=1.
        if (update->ntimestep == nevery){
            int idmc = maxcurvid(w);
            if(printflag==1 && myrank==0){write << update->ntimestep << " " << idmc << " " << topotype << endl;}
            placetopo(idmc);
        }
        
        //The next times it is called
        if (update->ntimestep > nevery){
            //Probability of moving
            double pmv = pmove();
            
            if(pmv<=prob){
                int idmc = maxcurvid(w);
                if(printflag==1 && myrank==0){write << update->ntimestep << " " << idmc << " " << topotype << endl;}
                changetype();
                placetopo(idmc);
            }
        }
    }
    
    next_reneighbor = update->ntimestep;
}


/*  -----------------------  */
/*  Choose randomly atom id  */
/*  -----------------------  */
int FixTopo2::randomid()
{   
    //The total number of atoms in the system
    bigint natoms = atom->natoms;
          
    //Choose Randomly a bead [1, natoms]
    int rng;
    rng = 1 + (rand() % natoms);
  
    return rng;
}


/*  -----------------------------------------------------*/
/*  Find the bead with maximum density (if flagtopo==3)  */
/*  -----------------------------------------------------*/
int FixTopo2::maxdensid(int w)
{
    /*  Gather data from all processors into one global array  */
    /***********************************************************/
   
    //We will use MPI_Gatherv for this purpose, so here is a description:
    //int MPI_Gatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, const int *recvcounts, const int *displs, MPI_Datatype recvtype, int root, MPI_Comm comm)

    //sendbuf   = **x:        starting address of send buffer (choice)
    //sendcount = nlocal:     number of elements in send buffer (integer)
    //sendtype  = MPI_DOUBLE: data type of send buffer elements (handle)
        
    //recvbuf    = xpos:             address of receive buffer (choice, significant only at root)
    //recvcounts = nlarray[mysize]:  integer array (of length group size) containing the number of elements that are received from each process (significant only at root)
    //displs     = disp[mysize]:     integer array (of length group size). Entry i specifies the displacement relative to recvbuf at which to place the incoming data from process i 
    //recvtype   = MPI_DOUBLE:       data type of recv buffer elements (significant only at root) (handle)
    //root       = myrank==0:        rank of receiving process (integer)
    //comm:                          communicator (handle)

    //Declaration of arrays used in this section
    int *nlarray = nullptr;
    int *disp    = nullptr;
    int *indice  = nullptr;
    double *xpos = nullptr;
    double *ypos = nullptr;
    double *zpos = nullptr;
    double **globalpos = nullptr;
  
    //Number of ranks used
    int mysize;
    MPI_Comm_size(MPI_COMM_WORLD, &mysize);
    
    //Rankid
    int myrank;
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    
    bigint natoms = atom->natoms;
    int nlocal = atom->nlocal;
    int nghost = atom->nghost;
    int m = nlocal+nghost;
    
    //Define in rank0 an array which entries are the number (nlocal), of atoms in each processor    
    if (myrank == 0) {
        memory->destroy(nlarray);
        memory->create(nlarray,mysize, "FixTopo2::maxdensid()");
    }
        
    MPI_Gather(&nlocal, 1, MPI_INT, nlarray, 1, MPI_INT, 0, MPI_COMM_WORLD);
        
    /* //Check   
    fprintf(screen,"mysize=%d, myrank=%d, nlocal=%d \n",mysize,myrank,nlocal);  
    if (myrank == 0) {
        for(int i=0; i<mysize; i++){
             fprintf(screen,"i=%d, nlarray[i]=%d \n",i, nlarray[i]);
        }
    }
    */
             
    int totlen = 0;
    //Space between two consecutive entries
    int space =0;
    //Define in rank0 an array which contains the displacement between entries, when using MPI_Gatherv.
    if (myrank == 0) {
        memory->destroy(disp);
        memory->create(disp,mysize, "FixTopo2::maxdensid()");
            
        disp[0]=0;
        totlen += nlarray[0];
        for (int i=1; i<mysize; i++) {
            totlen += nlarray[i];
            disp[i] = disp[i-1] + nlarray[i-1] + space;
        }
            
        // allocate memory for atom with positions
        memory->destroy(xpos);
        memory->destroy(ypos);
        memory->destroy(zpos);
        memory->destroy(indice);
        
        memory->create(xpos,natoms,"FixTopo2::maxdensid()");
        memory->create(ypos,natoms,"FixTopo2::maxdensid()");
        memory->create(zpos,natoms,"FixTopo2::maxdensid()");
        memory->create(indice,natoms,"FixTopo2::maxdensid()");
    }

    //We are going to gather the global indices that come from local in each rank
    tagint *tag = atom->tag;
    
    //Also we are going to gather the x,y,z position from each rank    
    double **x = atom->x;       
    double xtemp[nlocal];
    double ytemp[nlocal];
    double ztemp[nlocal];
    
    for(int l=0; l<nlocal; l++)
    {
        xtemp[l]=x[l][0];
        ytemp[l]=x[l][1];
        ztemp[l]=x[l][2];
    }
                
    MPI_Gatherv(tag, nlocal, MPI_INT, indice, nlarray, disp, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Gatherv(xtemp, nlocal, MPI_DOUBLE, xpos, nlarray, disp, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Gatherv(ytemp, nlocal, MPI_DOUBLE, ypos, nlarray, disp, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Gatherv(ztemp, nlocal, MPI_DOUBLE, zpos, nlarray, disp, MPI_DOUBLE, 0, MPI_COMM_WORLD);
                                   
    if (myrank == 0) {
        //Initialize array to store all the atom positions:
        memory->destroy(globalpos);
        memory->create(globalpos,natoms,3, "FixTopo2::maxdensid()");

        for(int i=0; i<natoms; i++){
            int g = indice[i];
            
            globalpos[g-1][0] = xpos[i];
            globalpos[g-1][1] = ypos[i];
            globalpos[g-1][2] = zpos[i];
            
            //fprintf(screen,"i=%d, x=%f, y=%f, z=%f \n", indice[i], xpos[i], ypos[i], zpos[i]);
        }
        
        //for(int i=0; i<natoms; i++){
        //    fprintf(screen,"i=%d, x=%f, y=%f, z=%f \n", i+1, globalpos[i][0], globalpos[i][1], globalpos[i][2]);
        //}
    }
    
    /* Local density is computed within a sphere of radius R */
    /*********************************************************/
    int *idmdarray = nullptr;
    int idmaxdens;

    if (myrank == 0) {
        memory->destroy(idmdarray);
        memory->create(idmdarray,mysize,"FixTopo2::maxdensid()");

        double R=pow(w,0.588);
        double Rsqd=R*R;
        int maxdens=0;
        int md;

        for(int g=0; g<natoms;g++){
            int localdens=0;

            double xtmp = globalpos[g][0];
            double ytmp = globalpos[g][1];
            double ztmp = globalpos[g][2];

            for (int j=0;j<natoms;j++){
                // calculate distance-squared between g,j atom-types
                double delx = xtmp - globalpos[j][0];
                double dely = ytmp - globalpos[j][1];
                double delz = ztmp - globalpos[j][2];
                double rsq = delx*delx + dely*dely + delz*delz;

                if (rsq < Rsqd) { localdens += 1;}
            }

            if(localdens>maxdens){
                maxdens=fabs(localdens);
                md=g+1;
            }
        }

        //We store a copy of the particle id with maximum density into this array
        //that is only defined in rank0
        for(int s=0; s<mysize; s++){
            idmdarray[s] = md;
        }
    }

    //Send the id with maximum density to all ranks
    MPI_Scatter(idmdarray, 1, MPI_INT, &idmaxdens, 1, MPI_INT, 0, MPI_COMM_WORLD);
    
    return idmaxdens;
}




/*  -----------------------------------------------------*/
/*  Find the bead with maximum curvature (if flagtopo==4)  */
/*  -----------------------------------------------------*/
int FixTopo2::maxcurvid(int w)
{
    /*  Gather data from all processors into one global array  */
    /***********************************************************/
   
    //We will use MPI_Gatherv for this purpose, so here is a description:
    //int MPI_Gatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, const int *recvcounts, const int *displs, MPI_Datatype recvtype, int root, MPI_Comm comm)

    //sendbuf   = **x:        starting address of send buffer (choice)
    //sendcount = nlocal:     number of elements in send buffer (integer)
    //sendtype  = MPI_DOUBLE: data type of send buffer elements (handle)
        
    //recvbuf    = xpos:             address of receive buffer (choice, significant only at root)
    //recvcounts = nlarray[mysize]:  integer array (of length group size) containing the number of elements that are received from each process (significant only at root)
    //displs     = disp[mysize]:     integer array (of length group size). Entry i specifies the displacement relative to recvbuf at which to place the incoming data from process i 
    //recvtype   = MPI_DOUBLE:       data type of recv buffer elements (significant only at root) (handle)
    //root       = myrank==0:        rank of receiving process (integer)
    //comm:                          communicator (handle)

    //Declaration of arrays used in this section
    int *nlarray = nullptr;
    int *disp    = nullptr;
    int *indice  = nullptr;
    double *xpos = nullptr;
    double *ypos = nullptr;
    double *zpos = nullptr;
    double **globalpos = nullptr;
  
    //Number of ranks used
    int mysize;
    MPI_Comm_size(MPI_COMM_WORLD, &mysize);
    
    //Rankid
    int myrank;
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    
    bigint natoms = atom->natoms;
    int nlocal = atom->nlocal;
    int nghost = atom->nghost;
    int m = nlocal+nghost;
    
    //Define in rank0 an array which entries are the number (nlocal), of atoms in each processor    
    if (myrank == 0) {
        memory->destroy(nlarray);
        memory->create(nlarray,mysize, "FixTopo2::maxcurvid()");
    }
        
    MPI_Gather(&nlocal, 1, MPI_INT, nlarray, 1, MPI_INT, 0, MPI_COMM_WORLD);
        
    /* //Check   
    fprintf(screen,"mysize=%d, myrank=%d, nlocal=%d \n",mysize,myrank,nlocal);  
    if (myrank == 0) {
        for(int i=0; i<mysize; i++){
             fprintf(screen,"i=%d, nlarray[i]=%d \n",i, nlarray[i]);
        }
    }
    */
             
    int totlen = 0;
    //Space between two consecutive entries
    int space =0;
    //Define in rank0 an array which contains the displacement between entries, when using MPI_Gatherv.
    if (myrank == 0) {
        memory->destroy(disp);
        memory->create(disp,mysize, "FixTopo2::maxcurvid()");
            
        disp[0]=0;
        totlen += nlarray[0];
        for (int i=1; i<mysize; i++) {
            totlen += nlarray[i];
            disp[i] = disp[i-1] + nlarray[i-1] + space;
        }
            
        // allocate memory for atom with positions
        memory->destroy(xpos);
        memory->destroy(ypos);
        memory->destroy(zpos);
        memory->destroy(indice);
        
        memory->create(xpos,natoms,"FixTopo2::maxcurvid()");
        memory->create(ypos,natoms,"FixTopo2::maxcurvid()");
        memory->create(zpos,natoms,"FixTopo2::maxcurvid()");
        memory->create(indice,natoms,"FixTopo2::maxcurvid()");
    }

    //We are going to gather the global indices that come from local in each rank
    tagint *tag = atom->tag;
    
    //Also we are going to gather the x,y,z position from each rank the unwrapped coordinates   
    double **x = atom->x;       
    double xtemp[nlocal];
    double ytemp[nlocal];
    double ztemp[nlocal];
    
    imageint *image = atom->image;
    double unwrap[3];
    
    for(int l=0; l<nlocal; l++)
    {
        domain->unmap(x[l],image[l],unwrap);
        
        xtemp[l]=unwrap[0];
        ytemp[l]=unwrap[1];
        ztemp[l]=unwrap[2];
    }
                
    MPI_Gatherv(tag, nlocal, MPI_INT, indice, nlarray, disp, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Gatherv(xtemp, nlocal, MPI_DOUBLE, xpos, nlarray, disp, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Gatherv(ytemp, nlocal, MPI_DOUBLE, ypos, nlarray, disp, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Gatherv(ztemp, nlocal, MPI_DOUBLE, zpos, nlarray, disp, MPI_DOUBLE, 0, MPI_COMM_WORLD);
                                   
    if (myrank == 0) {
        //Initialize array to store all the atom positions:
        memory->destroy(globalpos);
        memory->create(globalpos,natoms,3, "FixTopo2::maxcurvid()");

        for(int i=0; i<natoms; i++){
            int g = indice[i];
            
            globalpos[g-1][0] = xpos[i];
            globalpos[g-1][1] = ypos[i];
            globalpos[g-1][2] = zpos[i];
            
            //fprintf(screen,"i=%d, x=%f, y=%f, z=%f \n", indice[i], xpos[i], ypos[i], zpos[i]);
        }
        
        //for(int i=0; i<natoms; i++){
        //    fprintf(screen,"i=%d, x=%f, y=%f, z=%f \n", i+1, globalpos[i][0], globalpos[i][1], globalpos[i][2]);
        //}
    }


    /* Local curvature computed within a window */
    /********************************************/
    int *idmcarray = nullptr;
    int idmaxcurv;

    if (myrank == 0) {
        memory->destroy(idmcarray);
        memory->create(idmcarray,mysize,"FixTopo2::maxcurvid()");

        double maxcurv=0;
        double localcurv[natoms];
        int mc;

        for(int g=0; g<natoms;g++){
            localcurv[g]=0;

            //loop  over window centered at g
           	for (int j=0; j<=w;j++) {
               	double tang1[3];
               	double tang2[3];

                //Three consecutive beads
               	int b1=int(g+j-w/2.)-1;
           	    int b2=int(g+j-w/2.);
               	int b3=int(g+j-w/2.)+1;
           	    if(b1<0)b1=natoms+b1; // [0:natoms-1]
               	if(b2<0)b2=natoms+b2; // [0:natoms-1]
           	    if(b3<0)b3=natoms+b3; // [0:natoms-1]
               	if(b1>=natoms)b1=b1%natoms; // [0:natoms-1]
               	if(b2>=natoms)b2=b2%natoms; // [0:natoms-1]
           	    if(b3>=natoms)b3=b3%natoms; // [0:natoms-1]

                //Two consecutive tangents
               	for(int d=0; d<3; d++){
                   	tang1[d] = globalpos[b2][d]-globalpos[b1][d];
                   	tang2[d] = globalpos[b3][d]-globalpos[b2][d];
           	    }

               	double norm1=sqrt(tang1[0]*tang1[0]+tang1[1]*tang1[1]+tang1[2]*tang1[2]);
               	double norm2=sqrt(tang2[0]*tang2[0]+tang2[1]*tang2[1]+tang2[2]*tang2[2]);

               	double costheta=(tang1[0]*tang2[0]+tang1[1]*tang2[1]+tang1[2]*tang2[2])/(norm1*norm2);

           	    //integrate over window
               	//divide by w to get the "average"
               	localcurv[g]+=(1.0-costheta)*1.0/w;
            }

            if(localcurv[g]>maxcurv){
                maxcurv=localcurv[g];
                mc=g+1;
            }
        }

        //We store a copy of the particle id with maximum curvature into this array
        //that is only defined in rank0
        for(int s=0; s<mysize; s++){
            idmcarray[s] = mc;
        }
    }

    //Send the id with maximum density to all ranks
    MPI_Scatter(idmcarray, 1, MPI_INT, &idmaxcurv, 1, MPI_INT, 0, MPI_COMM_WORLD);

    return idmaxcurv;
}



/*  ------------------  */
/*  Moving probability  */
/*  ------------------  */
double FixTopo2::pmove()
{   
    //This is a number between 0 and 1.
    double pmv=rand()*1.0/RAND_MAX;
     
    return pmv;
}

/*  ----------  */
/*  Place topo  */
/*  ----------  */
void FixTopo2::placetopo(int id)
{
    //The total number of atoms in the system
    bigint natoms = atom->natoms;
        
    //The number  of owned atoms belonging to the core that executes this line
    int nlocal = atom->nlocal;

    //The number  of ghost atoms
    int nghost = atom->nghost;

    //Both owned and ghost atoms
    int n = nlocal+nghost;
    delete [] select;
    select = new int[n];

    //useful for converting from local to global
    tagint *tag = atom->tag;
  
    //Define the global id of the ltopo-consecutive beads that are going to change type.
    delete [] beads;
    beads = new int[ltopo];
    
    for (int i=0;i<ltopo;i++){
        int g=(id-ltopo+i);
        if(g<0)g+=natoms;
        beads[i]=g+1;
        //fprintf(screen,"Global id that will change type: %d \n",beads[i]);
    }

    //Firts select the local atoms that are going to change type according to the previously global selected ids
    for (int l = 0; l < nlocal; l++){
        select[l] = 0;
        for (int i=0; i<ltopo; i++){
            int g=beads[i];
            if (tag[l] == g) {
                select[l] = 1; 
                //fprintf(screen," l=%d, tag[l]=%d, g=%d \n",l,tag[l], g);
            }
        }
    }   
  
    //Change the type of previously selected atoms
    for (int l = 0; l < nlocal; l++) {
        if (!select[l]) continue;
        atom->type[l] = topotype;
    }
}



/*  ------------------------  */
/*  Change type intermediate  */
/*  ------------------------  */
void FixTopo2::changetype()
{    
    //The number  of owned atoms belonging to the core that executes this line
    int nlocal = atom->nlocal;

    //The number  of ghost atoms
    int nghost = atom->nghost;

    //Both owned and ghost atoms
    int n = nlocal+nghost;
    delete [] select;
    select = new int[n];

    //useful for converting from local to global
    tagint *tag = atom->tag;
    
    //Change the type of atoms from (topotype+1) to 1
    for (int l = 0; l < nlocal; l++) {
        if(atom->type[l] == topotype+1) atom->type[l] = 1;
    }
        
    //Change the type of atoms from (topotype) to (intertype)
    for (int l = 0; l < nlocal; l++) {
        if(atom->type[l] == topotype) atom->type[l] = intertype;
    }
}



/***********************************************************************/
/* Needed to write a restart file that can continue with the simulation*/
/***********************************************************************/
void FixTopo2::write_restart(FILE *fp)
{
  int n = 0;
  double list[4];
  list[n++] = random_equal->state();
  list[n++] = random_unequal->state();
  list[n++] = ubuf(next_reneighbor).d;
  list[n++] = ubuf(update->ntimestep).d;

  if (comm->me == 0) {
    int size = n * sizeof(double);
    fwrite(&size,sizeof(int),1,fp);
    fwrite(list,sizeof(double),n,fp);
  }
}

/* ----------------------------------------------------------------------
   use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTopo2::restart(char *buf)
{
  int n = 0;
  double *list = (double *) buf;

  seed = static_cast<int> (list[n++]);
  random_equal->reset(seed);

  seed = static_cast<int> (list[n++]);
  random_unequal->reset(seed);

  next_reneighbor = (bigint) ubuf(list[n++]).i;

  bigint ntimestep_restart = (bigint) ubuf(list[n++]).i;
  if (ntimestep_restart != update->ntimestep)
    error->all(FLERR,"Must not reset timestep when restarting fix topo2");
}