Commit2 JT 040119

Run this example as:

Run those examples as:

mpirun -np 3 lmp_g++ -partition 3x1 -in in.tad

mpirun -np 4 lmp_mpi -in in.gneb.iron -partition 4x1

You should be able to use any number of replicas >= 3.

Interpolate.x tries to perform a cubic polynomial interpolation

of the MEP found 

Compile the program with:

gcc interpolate.c -o interpolate.x -lm -lgsl

/* ------------------------------------------------------------------------

 Provide some explanation here

------------------------------------------------------------------------- */

/* ------------------------------------------------------------------------

   This program is a courtesy of Aleksei Ivanov (Univ. of Iceland)

   Contributing authors: Aleksei Ivanov (Univ. of Iceland), 

   			 Julien Tranchida (SNL)

------------------------------------------------------------------------- */

#include <stdio.h>

#include <math.h>

#include <time.h>

#include <string.h>

#include <stdlib.h>

// calculate cubic coefficients

void count_coefficient(double *V, double *F, double *R, double *a, double *b,

    double *c,double *d,int M){

  /* R = square of distance between images*/

  /* V = energy of images */

  /* F = projection of real real forces along the path? */

  int i;

  for(i = 0; i < M ; i++) {

    a[i] = ( -2.0*(V[i+1]-V[i])/R[i] - F[i] - F[i+1] ) / (R[i]*R[i]);

    //a[i] = ( -2.0*(V[i+1]-V[i])/sqrt(R[i]) - F[i] - F[i+1] ) / R[i];

    b[i] = ( 3.0*(V[i+1]-V[i]) + (2.0*F[i]+F[i+1])*R[i] ) / (R[i]*R[i]);

    //b[i] = ( 3.0*(V[i+1]-V[i]) + (2.0*F[i]+F[i+1])*sqrt(R[i]) ) / R[i];

    c[i] = -F[i];

    d[i] = V[i];

  }

}

// cubic spline

double spline(double a,double b,double c,double d,double x) {

  return (a*x*x*x + b*x*x + c*x + d);

}

int main() {

  int M=0; 			// M+1 = number of images

  double *fmdottan;		// projection of real forces on tangent path

  double *coords;		// initial value of reaction coords

  double *V; 			// energy of images

  double *dist;			// square of the distance between images

  double *a, *b, *c, *d ;	// coefficients of cubic functions 

  double x;			// reaction coordinate

  double A,B; 			// additional variables for rnd

  double length = 0.0;

  int i,j;

  FILE *data;

  printf("Enter M = number of images - 1 \n");

  scanf("%d",&M);

  // allocating and initializing tables

  a = calloc(M,sizeof(double));			// cubic coefficients

  b = calloc(M,sizeof(double));

  c = calloc(M,sizeof(double));

  d = calloc(M,sizeof(double));

  V = calloc((M+1),sizeof(double)); 		// energies

  coords = calloc((M+1),sizeof(double));	// reaction coordinates 

  fmdottan = calloc((M+1),sizeof(double));	// fm dot tangent

  dist = calloc(M+1,sizeof(double));		// distance between images

  // reading input file

  if((data=fopen("reac_coords_iron_verbose.dat","r")) == NULL) {

  //if((data=fopen("neb_init.dat","r")) == NULL) {

    printf("Incorrect input file name.");

    return 0;

  }

  for(j=0; j < M+1; j++) {

    fscanf(data,"%lf\t%lf\t%lf\t%lf\n",&coords[j],&V[j],&fmdottan[j],&dist[j]);

    length += dist[j];

    printf("%lf %lf %lf %lf\n",coords[j],V[j],fmdottan[j],dist[j]);

  }

  if( (fclose(data)) == 0) {

  	printf("Data stored, input file closed.\n ");

  }

  // calculate value of coefficients

  count_coefficient(V,fmdottan,dist,a,b,c,d,M);

  // plot result of the interpolation

  if( ( data=fopen("interpolation_result.dat","w") )== NULL) {

    printf("Interpolation file could not be open.");

    return 0;

  }

  A = B = 0.0;

  for(i = 0; i < M ; i++) {

    B += dist[i];

    printf("%13le\n",B);

    for(j = 0; j <= 1000; j++) {

      x = dist[i]*1.0e-3*j;

      fprintf(data,"%13lf\t%13le\n",(x+A)/length,spline(a[i],b[i],c[i],d[i],x));

    }

    A += dist[i];

  }

  return 0;

}

0.0000000	-3.915271	3.4995081e-17	2.4573077

0.14285714	-3.9148148	-0.00059075739	2.4573077

0.28571429	-3.9136926	-0.00072315767	2.4573077

0.42857143	-3.9127883	-0.0003191228	2.4573265

0.57142857	-3.9127883	0.0003191228	2.4403341

0.71428571	-3.9136926	0.00072315767	2.4044093

0.85714286	-3.9148148	0.00059075739	2.3766041

1.0000000	-3.915271	-4.1231828e-17	0.0000000