Bioinformatics tutorials: June 2013

Exported from Notepad++
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// input parameters

int np = 500;              // number of Ar particles
double density = 0.75;     // density of system
double temp = 1.0;         // desired temperature
double dt = 0.001;         // time step
double numstep = 1000;      // number of steps for MD

// constants and variables

double r[500][3];      // positions
double v[500][3];      // velocities
double f[500][3];          // forces

double L, t, temperature;  // box length
double ke, pe, total;      // kinetic energy, potential energy, total energy
int nc;                    // number of cells in each direction (cubic box)

// functions prototype

void initialize();
void calcforce();
void update();

int main () 
{
    int p, i;
    FILE *fh;
    fh = fopen("sim-box.xyz", "w");

    initialize();
    // writing initial configuration to file
    fprintf(fh, "%d\n\n", np);
    for (i = 0; i < np; i++) {
    fprintf(fh, "%s\t%f\t%f\t%f\n","Ar", r[i][0], r[i][1], r[i][2]);
    }

    calcforce();

    t = 0;
    for (t =0; t < numstep; t++) {
    update();
        // writing configuration at each step to file
        for (i = 0; i < np; i++) {
         fprintf(fh, "%s\t%f\t%f\t%f\n","Ar", r[i][0], r[i][1], r[i][2]);
        }
        ke = 0;
        for (p = 0; p < np; p++) {
        for (i = 0; i < 3; i++) {
            ke += 0.5 * v[p][i] * v[p][i];
        }
        }
        total = ke + pe;    
        printf("pe = %f ke = %f Total Energy = %f\n",pe, ke, total);
    }

    fclose(fh);
}

void initialize () {
    int i, j, k, m, n, p, nc, IX, IY, IZ, IREF, M;
    double L, CELL, CELL2, sumv2, sumv3, df, fs;
    double Sum[3]  = {0.0, 0.0, 0.0};

    // compute length of box
    L = powf( np / density, 1.0 / 3.0 );
    nc = powf( np / 4, 1.0 / 3.0 );

    /* CALCULATE THE SIDE OF THE UNIT CELL */
    CELL  = L / nc;
    CELL2 = 0.5 * CELL;
  
    /* BUILD THE UNIT CELL */
  
    /* SUBLATTICE A */
    r[0][0] =  0.0;
    r[0][1] =  0.0;
    r[0][2] =  0.0;
    /* SUBLATTICE B */
    r[1][0] =  0.0;
    r[1][1] =  CELL2;
    r[1][2] =  CELL2;
    /* SUBLATTICE C */
    r[2][0] =  CELL2;
    r[2][1] =  0.0;
    r[2][2] =  CELL2;
    /* SUBLATTICE D */
    r[3][0] =  CELL2;
    r[3][1] =  CELL2;
    r[3][2] =  0.0;
  
    /* CONSTRUCT THE LATTICE FROM THE UNIT CELL */
    M = 0;
    for( IZ=0; IZ<5; IZ++ )
        for( IY=0; IY<5; IY++ )
            for( IX=0; IX<5; IX++ ){
    
            for( IREF=0; IREF<4; IREF++ ){
      
                r[IREF+M][0] = r[IREF][0] + CELL * IX ;
                r[IREF+M][1] = r[IREF][1] + CELL * IY ;
                r[IREF+M][2] = r[IREF][2] + CELL * IZ ;
      
            }
        M = M + 4;
           }

    // random distributed initial velocities
    for (p = 0; p < np; p++)
    for (i = 0; i < 3; i++) {
        v[p][i] = (double)rand()/RAND_MAX - 0.5;
        Sum[i] += v[p][i];
    }

    // with zero total momentum
    for (p = 0; p < np; p++)
    for (i = 0; i < 3; i++)
        v[p][i] -= Sum[i] / np;

    // rescaled to desired temperature
    for (p = 0; p < np; p++)
    for (i = 0; i < 3; i++)
        sumv2 += v[p][i] * v[p][i];

    df = 3 * np - 3;
    fs = sqrt(temp * df / sumv2);

    printf("temperature before scaling= %f\n", sumv2/df);

    sumv3 = 0.0;
    for (p = 0; p < np; p++)
    for (i = 0; i < 3; i++) {
        v[p][i] *= fs;
            sumv3 += v[p][i] * v[p][i];
        }

    printf("temperature after scaling= %f\n", sumv3/df);

}

void calcforce () {
    int i, j, k, sign;
    double r2, rij[3], mag, L;
    double r2i, r6i, r12i, ff, fij[3];

    for (i = 0; i < np; i++)
    for (k = 0; k < 3; k++)
        f[i][k] = 0.0;

    pe = 0.0;
    L = powf(np / density, 1.0 / 3.0 );

    // loop over all pairs of particles
    for (i = 0; i < np - 1; i++) {
    for (j = i + 1; j < np; j++) {
        r2 = 0.0;
        for (k = 0; k < 3; k++) {
        rij[k] = r[i][k] - r[j][k];
        sign = rij[k] > 0.0 ? 1 : -1;
        mag = rij[k] * sign;
        if (mag > L / 2)         // use nearest image
            rij[k] -= L * sign;
        r2 += rij[k] * rij[k];
        }
        r2i = 1 / r2;
        r6i = r2i * r2i * r2i;
        r12i = r6i * r6i;
        ff = r2i * (48 * r12i - 24 * r6i);
        for (k = 0; k < 3; k++) {
        fij[k] = rij[k] * ff;
        f[i][k] += fij[k];
        f[j][k] -= fij[k];
        }
        pe += 4 * (r12i - r6i);
    }
    }
}

void update () {
    int p, k;

    // compute length of box
    L = powf( np / density, 1.0 / 3.0 );

    for (p = 0; p < np; p++)
    for (k = 0; k < 3; k++) {
        r[p][k] += v[p][k] * dt + 0.5 * f[p][k] * dt * dt;
        // impose periodic boundary conditions
        if (r[p][k] < 0.0)
        r[p][k] += L;
        if (r[p][k] >= L)
        r[p][k] -= L;
        v[p][k] += 0.5 * f[p][k] * dt;
    }
    calcforce();
    for (p = 0; p < np; p++)
    for (k = 0; k < 3; k++)
        v[p][k] += 0.5 * f[p][k] * dt;

}
Bioinformatics tutorials

Friday, June 28, 2013

Trajectory file conversion using wordom

Wednesday, June 12, 2013

MD simulation of LJ Particles (Ar) in C

About Me
