forked from lijiext/lammps
107 lines
3.2 KiB
C
107 lines
3.2 KiB
C
/* ------------------------------------------------------------------------
|
|
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 <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("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;
|
|
}
|