Created
September 29, 2019 04:27
-
-
Save anyuzx/fdd27b467c273665328955fbc111dfeb to your computer and use it in GitHub Desktop.
C++ code for Pivot algorithm https://www.guangshi.io/posts/pivot-algorithm/
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// | |
// This is a lattice random walk config generator using pivot algorithm. | |
// can be self-avoid or not self-avoid | |
// lattice_SAW(chain array, number of steps(monomers), step length(bond), self-avoid or no-self-avoid) | |
// | |
// Created by guangshi on 8/25/14. | |
// Copyright (c) 2014 Guang Shi. All rights reserved. | |
// | |
#include <math.h> | |
#include <random> | |
using namespace std; | |
void c_lattice_chain(double* chain, int N, double l0, int ve, int t){ | |
int i, j; | |
int acpt=0, duplicate; | |
int pick_pivot, pick_side, pick_rotate; | |
int symtry_oprtr[3][3]; | |
double c_x = 0, c_y = 0, c_z = 0; | |
// define the rotation matrix used in pivot algorithm | |
int rotate_matrix[9][3][3]={{{1,0,0},{0,0,-1},{0,1,0}},{{1,0,0},{0,-1,0},{0,0,-1}},{{1,0,0},{0,0,1},{0,-1,0}}, | |
{{0,0,1},{0,1,0},{-1,0,0}},{{-1,0,0},{0,1,0},{0,0,-1}},{{0,0,-1},{0,1,0},{-1,0,0}}, | |
{{0,-1,0},{1,0,0},{0,0,1}},{{-1,0,0},{0,-1,0},{0,0,1}},{{0,1,0},{-1,0,0},{0,0,1}} | |
}; | |
// use mt19937 random engine | |
std::random_device rd; | |
mt19937 re((unsigned int)rd()); | |
// set uniform integer distribution for choosing pivot point | |
std::uniform_int_distribution<int> ui_pivot(1,N-2); | |
// set uniform integer distribution for choosing side for rotating | |
std::uniform_int_distribution<int> ui_side(1,2); | |
// set unifrom interger distribution for choosing rotate matrix | |
std::uniform_int_distribution<int> ui_rotate(0,8); | |
// initialize the chain position | |
// Initially, chain is set to be straight along x-axis | |
for (i=0;i<N;i++){ | |
for (j=0;j<3;j++){ | |
if (j==0){ | |
chain[i*3+j] = double(i+1); | |
} | |
else{ | |
chain[i*3+j] = double(0); | |
} | |
} | |
} | |
// pivot algorithm. The number of successful rotating operations | |
// need to be the same as the number of steps | |
while (acpt <= t){ | |
// duplicate is used to check whether there are duplicate steps | |
// after rotating. | |
duplicate = 0; | |
// randomly pick the pivot point, side and rotate matrix | |
pick_pivot = ui_pivot(re); | |
pick_side = ui_side(re); | |
pick_rotate = ui_rotate(re); | |
// initialize the symmetry operator matrix | |
for (i=0;i<3;i++){ | |
for (j=0;j<3;j++){ | |
symtry_oprtr[i][j] = rotate_matrix[pick_rotate][i][j]; | |
} | |
} | |
// pick_side to rotate. | |
// pick_side == 1 means that the chain right to the pivot point is for rotating | |
if (pick_side == 1){ | |
// declare the old_chain and temp_chain | |
// old_chain is kept same. temp_chain is for rotating. | |
double old_chain [(pick_pivot+1)*3]; | |
double temp_chain [(N-pick_pivot-1)*3]; | |
// initilize old_chain and temp_chain | |
for (i=0;i<N;i++){ | |
for (j=0;j<3;j++){ | |
if (i<=pick_pivot){ | |
old_chain[i*3+j]=chain[i*3+j]; | |
} | |
else if(i>pick_pivot){ | |
temp_chain[(i-pick_pivot-1)*3+j]=chain[i*3+j]; | |
} | |
} | |
} | |
// apply the rotation to the chain | |
int new_chain[(N-pick_pivot-1)*3]; | |
for (i=0;i<N-pick_pivot-1;i++){ | |
for (j=0;j<3;j++){ | |
new_chain[i*3+j] = (temp_chain[i*3+0]-chain[pick_pivot*3+0])*symtry_oprtr[j][0]+(temp_chain[i*3+1]-chain[pick_pivot*3+1])*symtry_oprtr[j][1]+(temp_chain[i*3+2]-chain[pick_pivot*3+2])*symtry_oprtr[j][2]+chain[pick_pivot*3+j]; | |
} | |
} | |
// check whether it is self-avoiding or not | |
if (ve == 1){ | |
// looping the old_chain and new_chain to check duplicates | |
for (i=0;i<N-pick_pivot-1;i++){ | |
for (j=0;j<=pick_pivot;j++){ | |
if(new_chain[i*3+0]==old_chain[j*3+0]&&new_chain[i*3+1]==old_chain[j*3+1]&&new_chain[i*3+2]==old_chain[j*3+2]){ | |
duplicate = 1; | |
break; | |
} | |
} | |
if(duplicate == 1){break;} // there are duplicates, break the loop | |
} | |
if(duplicate == 0){ | |
for(i=pick_pivot+1;i<N;i++){ | |
for(j=0;j<3;j++){ | |
chain[i*3+j] = new_chain[(i-pick_pivot-1)*3+j]; | |
} | |
} | |
acpt++; | |
continue; // if there is no duplicate steps, then count the number of succesful rotating steps plus 1 | |
} | |
else if(duplicate == 1){ | |
continue; // if there are duplicate steps, then start the pivot algorithm from the beginning | |
} | |
} | |
// no self-avoiding | |
else{ | |
for(i=pick_pivot+1;i<N;i++){ | |
for(j=0;j<3;j++){ | |
chain[i*3+j] = new_chain[(i-pick_pivot-1)*3+j]; | |
} | |
} | |
acpt++; | |
continue; | |
} | |
} | |
// pick_side == 2 means that the chain left to the pivot point is for rotating | |
else if(pick_side == 2){ | |
double old_chain [(N-pick_pivot)*3]; | |
double temp_chain [(pick_pivot)*3]; | |
for (i=0;i<N;i++){ | |
for (j=0;j<3;j++){ | |
if (i<pick_pivot){ | |
temp_chain[i*3+j]=chain[i*3+j]; | |
} | |
else if(i>=pick_pivot){ | |
old_chain[(i-pick_pivot)*3+j]=chain[i*3+j]; | |
} | |
} | |
} | |
// apply the rotate to the chain | |
int new_chain[pick_pivot*3]; | |
for (i=0;i<pick_pivot;i++){ | |
for (j=0;j<3;j++){ | |
new_chain[i*3+j] = (temp_chain[i*3+0]-chain[pick_pivot*3+0])*symtry_oprtr[j][0]+(temp_chain[i*3+1]-chain[pick_pivot*3+1])*symtry_oprtr[j][1]+(temp_chain[i*3+2]-chain[pick_pivot*3+2])*symtry_oprtr[j][2]+chain[pick_pivot*3+j]; | |
} | |
} | |
if (ve == 1){ | |
for (i=0;i<pick_pivot;i++){ | |
for (j=0;j<N-pick_pivot;j++){ | |
if(new_chain[i*3+0]==old_chain[j*3+0]&&new_chain[i*3+1]==old_chain[j*3+1]&&new_chain[i*3+2]==old_chain[j*3+2]){ | |
duplicate = 1; | |
break; | |
} | |
} | |
if(duplicate == 1){break;} | |
} | |
if(duplicate == 0){ | |
for(i=0;i<pick_pivot;i++){ | |
for(j=0;j<3;j++){ | |
chain[i*3+j] = new_chain[i*3+j]; | |
} | |
} | |
acpt++; | |
continue; | |
} | |
else if(duplicate == 1){ | |
continue; | |
} | |
} | |
else{ | |
for(i=0;i<pick_pivot;i++){ | |
for(j=0;j<3;j++){ | |
chain[i*3+j] = new_chain[i*3+j]; | |
} | |
} | |
acpt++; | |
continue; | |
} | |
} | |
} | |
// calculate the center of the walk(chain) | |
for (i=0;i<N;i++){ | |
c_x = c_x + chain[i*3]; | |
c_y = c_y + chain[i*3+1]; | |
c_z = c_z + chain[i*3+2]; | |
} | |
c_x = c_x/N; | |
c_y = c_y/N; | |
c_z = c_z/N; | |
int c_r[3] = {int(c_x),int(c_y),int(c_z)}; | |
// loop over the chain. move the chain to center and mulitply by step size. | |
for (i=0;i<N;i++){ | |
for (j=0;j<3;j++){ | |
chain[i*3+j] = l0*(chain[i*3+j]-c_r[j]); | |
} | |
} | |
return; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment