main.cpp.dump

/*
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>

using namespace std;
# define PI           3.14159265358979323846  */
/* pi *//*


double  *UAvg_A;
double *Var_A;
const int MCSteps = 100;
long l = -1;
// simulation parameters

const int N = 32;               // number of particles
const double rho = 1.0;         // density (number per unit volume)
double lambda;            //order parameter computed from atom configuration
double L;                 // size of whole cube
double a;               //LatticeA constant
double M;                //number of cubes in LatticeA
double **LatticeA;     // positions,
double **LatticeB;     // positions,
double rCutOff;     // cut-off on Lennard-Jones potential

double ran0(long *);
void Initialize();

double MinImage(int, int, string flag);
void computeOrderParameter();

void MMSteps(int, double, string filename, double step, string flag);
double computeEnergy(string flag);
void writeHeaders(ofstream& f, int index, string *fnames);

int main() {
    time_t s_i;
    s_i = time(NULL);
    string filenames_A[3] = {"fcc_a0.txt", "fcc_a1.txt", "fcc_a2.txt"};
    string filenames_B[3] = {"fcc_b0.txt", "fcc_b1.txt", "fcc_b2.txt"};
    if (remove("fcc_a0.txt") == 0
        && remove("fcc_a1.txt") == 0
        && remove("fcc_a2.txt") == 0
        && remove("fcc_a2.txt") == 0
        && remove("fcc_a2.txt") == 0
        && remove("fcc_a2.txt") == 0
            )
        puts("Files successfully deleted");
    else perror("Error deleting file");
    double temperatures[3] = {0.1, 0.5, 1};
    double steps[3] = {0.06, 0.065, 0.08};
    ofstream file_a;
    ofstream file_b;

    for (int i = 0; i < 3; i++) {
        writeHeaders(file_a, i, filenames_A);
        writeHeaders(file_b, i, filenames_B);
        Initialize();
        for (int s = 0; s < MCSteps; s++) {
            MMSteps(s, temperatures[i], filenames_A[i], steps[i], "A");
            MMSteps(s, temperatures[i], filenames_A[i], steps[i], "B");
        }
    }

    delete [] UAvg_A;
    delete [] Var_A;
    delete [] LatticeA;
    delete [] LatticeB;

    time_t s_f;
    s_f = time(NULL);
    cout << "Runtime : " << s_f - s_i <<" seconds \n";
    system ("pause");
    return 0;
}

void computeOrderParameter(string flag){
    //compute the order parameter after each MonteCarlo time step

    lambda = 0;

    if (flag.compare("A") == 0) {
        for (int i = 0; i < N; i++){

            lambda += cos((4 * PI * LatticeA[i][0]) / a) +
                      cos((4 * PI * LatticeA[i][1]) / a) +
                      cos((4 * PI * LatticeA[i][2]) / a);
        }
    } else {
        for (int i = 0; i < N; i++){

            lambda += cos((4 * PI * LatticeB[i][0]) / a) +
                      cos((4 * PI * LatticeB[i][1]) / a) +
                      cos((4 * PI * LatticeB[i][2]) / a);
        }
    }
    lambda /= -3 * (double) N;
}

double MinImage(int i, int j, string flag) {
    // find separation using closest image convention
    double dr[3];
    double sum = 0;
    for (int d = 0; d < 3; d++) {
        if(flag.compare("A") == 0 ){
            dr[d] = LatticeA[i][d] - LatticeA[j][d];
        }
        else if(flag.compare("B") == 0){
            dr[d] = LatticeB[i][d] - LatticeB[j][d];
        }
        else{
            dr[d] = LatticeA[i][d] - LatticeA[j][d];
        }

        if (dr[d] >= 0.5*L) dr[d] -= L;
        if (dr[d] < -0.5*L) dr[d] += L;

        sum += dr[d] * dr[d];
    }
    return sum;
}

double computeEnergy(string flag) {
    double U = 0;
    for (int i = 0; i < N-1; i++)  {             // all distinct pairs
        for (int j = i+1; j < N; j++) {         // of particles i,j
            double rSqd = MinImage(i, j, flag);
            if (rSqd < rCutOff*rCutOff){
                U += (pow(1 / rSqd, 6) - pow(1 / rSqd, 3));
            }
        }
    }
    return 4 * U;
}

double *Move(int index, double dr, string flag) {
    static double delta[3] = {};
    double t = ran0(&l);
    t = 2 * t - 1;
    t = acos(t);
    double ph = ran0(&l);

    delta[0] = dr * cos(t) * sin(2 * ph * PI);
    delta[1] = dr * sin(t) * sin(2 * ph * PI);
    delta[2] = dr * cos(2 * ph * PI);

    if (flag.compare("A") == 0) {
        LatticeA[index][0] += delta[0];
        LatticeA[index][1] += delta[1];
        LatticeA[index][2] += delta[2];
    } else {
        LatticeB[index][0] += delta[0];
        LatticeB[index][1] += delta[1];
        LatticeB[index][2] += delta[2];
    }
    return delta;
}

void UnMove(int index, double *delta, string flag) {

    if (flag.compare("A") == 0) {
        LatticeA[index][0] -= delta[0];
        LatticeA[index][1] -= delta[1];
        LatticeA[index][2] -= delta[2];
    } else {
        LatticeB[index][0] -= delta[0];
        LatticeB[index][1] -= delta[1];
        LatticeB[index][2] -= delta[2];
    }
}
//function to apply the MonteCarlo Metropolis Step
void MMSteps(int s, double T, string filename, double step, string flag) {
    ofstream file;
    file.precision(8);
    file.open(filename, ios::app);

    long ll = -1;
    double Ui = 0;
    double U = 0;
    double delta_r = step; //0.03;
    double AcceptanceRatio = 0;
    computeOrderParameter();

    vector<int> indexes;
    for (int i = 1; i < N; ++i) indexes.push_back(i);

    while (indexes.size() != 0) {
        // choose a random point, here we are using idx as index
        random_shuffle(indexes.begin(), indexes.end());
        int idx = indexes[indexes.size() - 1];
        indexes.pop_back();
        double Uold = computeEnergy(flag);

        double *delta= Move(idx, delta_r, flag);
        // Apply periodic boundary conditions
        for (int k = 0; k < 3; k++){
            if (LatticeA[idx][k] < 0) LatticeA[idx][k] += L;
            if (LatticeA[idx][k] >= L)LatticeA[idx][k] -= L;
            if (LatticeB[idx][k] < 0) LatticeB[idx][k] += L;
            if (LatticeB[idx][k] >= L)LatticeB[idx][k] -= L;

        }

        //check energy difference due to move
        double Unew = computeEnergy(flag);
        double delta_U = Unew - Uold;
        if (delta_U > 0 && ran0(&ll) > exp(-(delta_U / T))) { //move is rejected
            UnMove(idx, delta,flag);
            Ui = Uold;
        }
        else {//move is accepted
            Ui = Unew;
            AcceptanceRatio +=1;
        }
        U += (Ui)/((double)N);
    }

    //compute the average energy and standard deviation using the method on
    UAvg_A[s] = s == 0 ? U : UAvg_A[s - 1] + ((U - UAvg_A[s - 1]) / s);

    Var_A[s] = s == 0 ? 0 : Var_A[s - 1] * s + (U - UAvg_A[s - 1]) * (Ui - UAvg_A[s]);
    Var_A[s] /= (s + 1);
    computeOrderParameter();
    AcceptanceRatio = (AcceptanceRatio * 100)/(double)N;

    file << s << "\t " << setw(12) << U << "\t            ";
    file << setw(12) << UAvg_A[s] << "\t  ";//running average
    file << setw(12) << Var_A[s] << "\t  ";//running variance
    file << setw(12) << lambda << "\t ";
    file << setw(12) << AcceptanceRatio << "\t";
    file << setw(12) << T << "\n";
    file.close();
}

void Initialize() {
    LatticeA = new double *[N];
    LatticeB = new double *[N];
    for (int i = 0; i < N; i++) {
        LatticeA[i] = new double[3];
        LatticeB[i] = new double[3];
    }

    UAvg_A = new double [MCSteps];
    Var_A = new double [MCSteps];

    // compute side of cube from number of particles and number density
    L = pow(N / rho, 1.0/3);
    rCutOff = 0.49 * L;
    // find M large enough to fit N atoms on an fcc LatticeA
    M = 1;
    while (4 * M * M * M < N)
        ++M;
    a = L / M;           // LatticeA constant of conventional cell
    double dx[4] = {0.0, 0.5, 0.5, 0.0};//
    double dy[4] = {0.0, 0.0, 0.5, 0.5};//
    double dz[4] = {0.0, 0.5, 0.0, 0.5};//

    int n = 0;                  // atoms placed so far
    for (int x = 0; x < M; x++){
        for (int y = 0; y < M; y++){
            for (int z = 0; z < M; z++){
                for (int p = 0; p < 4; p++){
                    if (n < N) {
                        LatticeA[n][0] = (0.25 + x + dx[p]) * a - 0.5*L;
                        LatticeA[n][1] = (0.25 + y + dy[p]) * a;
                        LatticeA[n][2] = (0.25 + z + dz[p]) * a - 0.5*L;

                        LatticeB[n][0] = (0.25 + x + dx[p]) * a - 0.5*L;
                        LatticeB[n][1] = n==0?
                                         -(0.75 + y + dy[p]) * a :
                                         -(0.25 + y + dy[p]) * a;
                        LatticeB[n][2] = (0.25 + z + dz[p]) * a - 0.5*L;
                        ++n;
                    }//end if
                }//end for p
            }//end for z
        }//end for y
    }//end for x
}

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876
double ran0(long *idum) {
    long     k;
    double   ans;
    *idum ^= MASK;
    k = (*idum)/IQ;
    *idum = IA*(*idum - k*IQ) - IR*k;
    if(*idum < 0) *idum += IM;
    ans=AM*(*idum);
    *idum ^= MASK;
    return ans;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef MASK
#undef PI

void writeHeaders(ofstream& f, int index, string *fnames){
    f.open(fnames[index], ios::trunc);
    f.precision(8);
    f << "MCStep s " << setw(12) << "\t " <<
    "Energy_U" << "\t" <<
    setw(12) << "Average_Energy_<U>" << "\t  " <<
    setw(12) << "Variance" << "\t " <<
    setw(12) << "Order_Parameter" << "\t " <<
    setw(12) << "Acceptance_Ratio" << "\t " <<
    setw(12) << "Temperature\n";
}*/