/* ----------------------------------------------------------------------
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");
}